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ABSTRACT 


This  document  includes  a report  that  describes  the  theoretical  basis  of  the  program 
DTAM1  and  a users  manual  for  the  program. 

DTAM1  is  a general  purpose  building  energy  simulation  program  that  was  developed  to 
demonstrate  an  approach  to  building  energy  simulation  based  upon  discrete  analysis 
techniques,  including,  but  not  limited  to,  the  Finite  Element  Method,  used  in  other  fields  of 
physical  simulation.  It  is  the  product  of  a first  phase  of  development  of  Discrete  Thermal 
Element  Analysis  Techniques  for  Building  Energy  Simulation  that  are  expected  to  provide  a 
means  to  unify  existing  building  energy  simulation  theory. 

DTAM1  provides  a library  of  discrete  thermal  elements,  that  may  be  assembled  to  model 
thermal  systems  idealized  to  have  constant  material  and  heat  transfer  properties  (i.e.,  linear 
idealizations),  including: 

- ID  two-node  thermal  resistance  elements 

- single-node  lumped  capacitance  elements 

- two-node  fluid  flow  loop  element 

- ID  two-to-four  node  isoparametric  conduction  Finite  Elements 

- 2D  four-node  isoparametric  conduction  Finite  Elements  (planar  and  axisymmetric) 

Equations  defining  a variable  node  mean  radiant  temperature  element  are  also  presented  in 
the  report. 

Steady  state  and  transient  analysis  capabilities  are  included.  Temperature,  heat  flow  rate, 
and  convective  boundary  conditions  may  be  modeled  and  system  temperature  variables  may 
be  constrained  to  be  equal  so  that  mixed  assemblages  of  ID  and  2D  elements  may  be 
employed. 


KEY  WORDS:  building  energy  simulation,  building  dynamics,  computer  simulation 
techniques,  discrete  analysis  techniques,  discrete  thermal  elements,  dynamic  simulation,  finite 
element  analysis,  DTAM1 
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REPORT 


1.  Introduction 

DTAM1  is  a general  purpose  building  energy  simulation  program  that  was  developed  to 
demonstrate  an  approach  to  building  energy  simulation  based  upon  discrete  analysis 
techniques,  including,  but  not  limited  to,  the  Finite  Element  Method,  used  in  other  fields  of 
physical  system  simulation.  It  is  the  product  of  a first  phase  of  development  of  Discrete 
Thermal  Element  Analysis  Techniques  for  Building  Energy  Simulation  that  are  expected  to 
provide  a means  to  unify  existing  building  energy  simulation  theory.  Portions  of  the 
program  are  based  on  the  Finite  Element  conduction  analysis  program,  *HEAT*  .written  by 
Robert  Taylor,  Department  of  Civil  Engineering,  U.  C.  , Berkeley  [1]  . DTAM1  was 
developed  using  the  CAL-SAP  library  of  subroutines  developed  by  Ed  Wilson,  Department 
of  Civil  Engineering,  U.C.,  Berkeley  [2]  . 

DTAM1  provides  a library  of  discrete  thermal  elements,  that  may  be  assembled  to  model 
thermal  systems  idealized  to  have  constant  material  and  heat  transfer  properties  (i.e.,  linear 
idealizations),  including: 

- ID  two-node  thermal  resistance  elements 

- single-node  lumped  capacitance  elements 

- two-node  fluid  flow  loop  elements 

- ID  two-to-four  node  isoparametric  conduction  Finite  Elements 

- 2D  four-node  isoparametric  conduction  Finite  Elements  (planar  and  axisymmetric) 

Equations  defining  a variable  node  mean  radiant  temperature  element  are  also  presented  in 
the  report. 

Static  and  dynamic  analysis  capabilities  are  included.  Temperature,  heat  flow  rate,  and 
convective  boundary  conditions  may  be  modeled  and  system  temperature  variables  may  be 
constrained  to  be  equal  so  that  mixed  assemblages  of  ID  and  2D  elements  may  be 
employed. 


2.  Theoretical  Basis 
2.1  Thermal  Model 

A building  thermal  system  may  be  modeled  as  an  assemblage  of  discrete  thermal  elements 
connected  at  discrete  system  nodes.  Typically,  these  nodes  will  be  associated  with  discrete 
points  in  the  physical  system,  although  this  is  not,  in  general,  necessary.  It  will  be  shown 
that  system  equations,  governing  the  behavior  of  the  building  thermal  system,  may  be 
assembled  from  element  equations  that  govern  the  behavior  of  individual  elements. 

A hypothetical  example  is  illustrated  below.  Here,  a simple  house  is  idealized  to  model 
dynamic  thermal  exchanges  between  the  house  and  its  environment  using  resistance 
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elements,  lumped  capacitance  elements,  and  ID  and  2D  conduction  finite  elements. 
Idealization  of  a thermal  system  as  an  assemblage  of  discrete  thermal  elements  is  facilitated 
by  the  use  of  diagrams  such  as  the  one  illustrated  below.  The  process  of  mathematical 
assembly  of  discrete  element  equations  is  analogous  to  the  graphical  assembly  of  element 
icons  used  in  such  diagrams  and  provides  the  analyst  with  an  intuitive  aide  to  the  process  of 
thermal  idealization. 


System  Node 


Lumped 

Capacitance 


ID  2-Node  _ 
Finite  Element 


ID  2-Node  Finite  Element 
MRT  Sub-Network 

Simple  Resistance  Elems. 
q j Temp.  DOF 


Constrained  to  be  Equal 
2D  4-Node  Finite  Element 


Fia.  1 Hypothetical  Example  of  a Thermal  Model  for  a Residence 


2.2  System  Degrees  of  Freedom 

With  each  system  node  we  associate  a temperature  degree  of  freedom  (DOF),  Tj , and  a 
direct  heat  flow  rate  DOF,  Q j (i.e.,  heat  flow  rate  directly  from  a source  into  the  system 
node).  These  variables  define  the  state  of  the  system,  at  any  time,  and  will  be  referred  to  as 
the  system  degrees  of  freedom  . The  system  DOFs  may  be  represented  as  vectors  for 
symbolic  mathematical  manipulation: 

{T)-{Ti,T2,T3 Tn}T  ; the  system  temperature  vector  (2.1) 

{Q}  = { Qi  , 02 , 03 Qn  }T  ; the  system  direct  heat  flow  rate  vector  (2.2) 

for  a system  with  n system  nodes.  The  system  temperature  vector  will  be  treated  as  the 
principal  dependent  variable. 

2.3  General  Form  of  Element  Equations 

One  or  more  nodes  may  be  associated  with  each  element.  With  each  element  node,  i,  we 

associate  an  element  temperature  DOF,  Tej,  and  net  heat  flow  rate  DOF,  qenet.j  (i.e.,  heat 
flow  rate  from  the  element  node  into  the  element): 
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These  element  DOFs  may  also  be  represented  as  vectors  for  symbolic  mathematical 
manipulation: 

{T6}  = { T®i  , T®2  > T®3  , , T®m  }T  ; the  element  temperature  vector  (2.3) 

{qenet}  = f qenet-1  ’ ^061-2 ^net-m  }T  : the  element  net  heat  flow  rate  vector  (2-4) 

for  an  element  with  m nodes. 

In  general,  we  attempt  to  describe  the  overall  energy  balance  of  an  element  by  equations 
of  the  form: 

{qenet}  - [ke](Te}  + [ce]{dTe/dt}  - {qe}  (2.5) 

where: 

[ke]  = matrix  of  element  conductance  transfer  coefficients 
= the  element  conductance  matrix 

= [ke(t,  Te)] , in  general 

[ce]  = matrix  of  element  thermal  capacity  coefficients;  the  element  capacitance 
matrix 

= [c®(t,T®)] , in  general 

{q®}  = vector  of  element-generated  heat  flow  rates;  the  element  internal  generation 

rate  vector 

For  elements  employed  in  the  program  DTAM1,  the  element  conductance  matrices,  [ke], 

will  be  positive  semidefinite  and  the  element  capacitance  matrices,  [c®],  will  be  positive 
definite. 

Physically,  the  terms  on  the  right  hand  side  of  equation  (2.5)  will  often  represent  heat 
flow  rate  into  the  element,  at  each  node,  by  conduction,  convection  and/or  radiation,  heat 
flow  into  storage,  and  heat  flow  generated  internally  within  the  element,  respectively. 
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Other  forms  of  element  equations  could  have  been  considered,  but  equations  of  this  form 
embrace  a large  variety  of  thermal  models  and  lead  to  systems  of  equations  that  have  been 
extensively  studied  and  for  which  many  numerical  solution  strategies  have  been  published. 
(Other  forms  based  upon  conduction  transfer  functions,  transmission  matrix  formulations, 
and  equations  associated  with  fluid  flow  heat  transfer  phenomena  have  been  formulated  and 
are  presently  under  study.) 


2.4  System  Equations 

Requiring  conservation  of  energy  at  each  system  node  we  demand: 


{ (direct  heat  flow  rate)  = X (net  heat  f|ow  rate  int0  element ) Jsystem 

connected  nocie 

elements 

(2.6) 

or,  for  an  arbitrary  system  node,  n,  with  connected  elements  "a",  "b",  "c",  ... 

Qn  - 1an.,-k  + 1bnet-l  + ^net-i  + --  <2'7> 

where  k,  I,  i,  ...  are  the  element  nodes  of  elements  a,  b,  c respectively,  that  are 

associated  with  the  system  node  n. 


Fig.  3 Conservation  of  Energy  at  System  Node 

If  individual  element  equations,  of  the  form  of  equation  (2.5),  are  substituted  directly  into 
the  equilibrium  relation  of  equation  (2.7)  a system  of  equations  that  relate  system  direct  heat 
flow  rate  DOFs,  (Q)  , to  the  collection  of  element  temperature  DOFs,  {T3},  {Tb},  {TC},  ... 
would  result.  As  each  element  temperature  DOF  is  associated  with  a specific  system 
temperature  DOF  it  is  possible  to  reform  these  equations  into  the  so-called  system 
equations  , that  have  the  form: 

[ K ]{  T } + [ C ]{dT  /dt)  = { E } (2.8) 
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where; 

[ K ] = matrix  of  system  conductance  transfer  coefficients;  the  system  conductance 

matrix 

[ C ] = matrix  of  system  thermal  capacity  coefficients;  the  system  capacitance  matrix 

{ E } = vector  of  source/sink  flow  rates;  the  system  excitation  vector 

For  elements  employed  in  the  program  DTAM1,  the  system  conductance  matrix,  [K],  will  be 
positive  semidefinite  and  the  system  capacitance  matrix,  [C],  will  be  positive  definite. 

The  process  of  forming  the  system  equations  from  the  element  equations  is  known  as  the 
element  assembly  process  . It  may  be  represented  formally  by  first  establishing  the  one-to- 
one  correspondence  between  each  element’s  DOFs  and  the  system  DOFs  through  simple 
Boolean  Transformations  of  the  form: 

(T*)  = [ Be]  {T)  (2.9a) 

or,  equivalently,  as  the  correspondence  is  one-to-one: 

{T}  = [ B®]T  {T«}  (2.9b) 

where: 

[ Be]  = the  transformation  matrix  for  element  "e"  ; a Boolean  transformation  matrix 
consisting  of  zeros  and  ones  as  the  element  DOFs  are  either  equal  to  a system 
DOF  (i.e.,  1)  or  not  (i.e.,  0). 

The  same  correspondence  exists  between  each  element’s  net  heat  flow  rate  DOFs  and  the 
contribution  of  that  element  to  the  equilibrium  at  each  of  the  system  nodes,  represented, 

here,  by  {Qe}  , as: 

(Qe)  = [ Be]T  {Qenet}  <2.10) 

Rewriting  equation  (2.7)  in  terms  of  the  element  contributions,  (Qe) , for  all  system  nodes: 

{0}  = {Qa}  + {Qb}  + {Qc}  + ...  (2.11a) 

or 

{Q}  = I{Q6}  (2.11b) 

a,  b,  c,  ... 

and  substituting  the  transformation  expressions  (2.9a),  (2.10),  and  (2.5)  we  obtain 
equation  (2.8),  above,  with: 
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(2.12) 


[ K ] = [Ba]T[ka][Ba]  + [Bb]T[kb][Bb]  + [Bc]T[kc][Bc]  + ... 

[ C ] = [Ba]T[ca][Ba]  + [Bb]T[cb][Bb]  + [B^T^KB0]  + ...  (2.13) 

{ E } . {Q}  + { [Ba]T{qa}  + [B^qb)  + [B^qCj  + ...  (2.14) 

The  formal  representation  of  the  assembly  process,  equations  (2.12)  , (2.13)  , and  (2.14), 
is  useful  for  theoretical  consideration.  Practically,  however,  the  system  matrices  may  be 
assembled  more  directly  and  efficiently  using  the  so-called  "LM  Algorithm"  [3]  . This 
algorithm  is  used  to  assemble  element  equations  to  form  the  system  equations  in  the 
program  DTAM1. 


2.5  Boundary  Conditions 


Some  nodal  temperatures,  {Tp},  may  be  prescribed  (e.g.,  outside  air  temperature),  while 
the  rest,  (T{)  may  be  considered  variable  or  free.  At  those  nodes  where  nodal  temperatures 

are  variable,  direct  heat  flow  rate  (or  in  the  case  of  some  of  the  elements,  convective  flow 
rates)  may  be  prescribed.  Equation  (2.8)  , then  , may  be  partitioned  as: 


i 

Kfp' 

( ) + 

i 

o 

1 

CL 

o 

Q. 

* 

1 

Kpp 

lTp| 

1 

o 

T> 

1 

CL 

CL 

O 

dll) 


dt  / 


(Note:  for  lumped  mass  idealizations  [Cfp]  = [Cpf]  = [0] .) 


(2.15) 


The  first  equation  of  equations  (2.15)  provides  the  governing  equation  of  the  free 
response  of  the  system: 

[KffHTf}  + [C„]{dT f/dt}  = (Ef)  - [K{p](Tp}  - [Cfp]{dTp/dt}  (2.16) 


The  right  hand  side  of  this  equation  defines  an  effective  excitation,  that  includes  the  effect 
of  prescribed  temperatures  on  system  response. 

Equation  (2.16)  may  be  solved  by  a variety  of  methods  and  the  solution  for  (Tf)  may  then 
be  substituted  into  the  second  equation  above: 

{Ep}  = [Kpf]{Tf)  + [KppHTp}  + [CpfHdTf/dt}  + [CppKdTp/dt)  (2.17) 


to  determine  the  excitation  quantities,  {Ep},  needed  to  maintain  the  prescribed  temperatures, 
{Tp} . 

A prescribed  nodal  temperature,  Tj  , may  also  be  imposed  numerically  by  scaling  the 
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diagonal  terms  corresponding  to  the  prescribed  temperature  DOF  of  either  the  system 
conductance  matrix  - term  Kjj  - and/or  the  system  capacitance  matrix  - term  Cjj  - by  a large 
value  and  setting  the  corresponding  excitation  term,  Ej  , to  the  product  of  the  prescribed 
temperature  and  the  appropriate  combination  of  Kjj  and  Cjj  (depending  upon  the  solution 
type  employed,  steady,  harmonic,  or  predictor-corrector  integration)  . This  method  is  simple, 
effective,  and  easily  implemented.  It  is,  therefore,  employed  in  DTAM1. 


2.6  Solution  of  System  Equations 

The  solution  of  equation  (2.8)  for  (T(t)}  defines  the  system  response  to  a given 
excitation,  (E(t)}  Response  analysis  may  be  classified  by  the  nature  of  the  excitation,  the 
nature  of  the  system  (i.e.,  linear  or  nonlinear;  constant  property  or  nonconstant  property), 
and  the  type  and  characteristics  of  the  analytical  or  numerical  method  used  to  obtain  the 
solution.  Here  we  shall  limit  consideration  to  linear  systems  with  constant  properties, 
subjected  to  either  steady  excitation,  harmonic  excitation,  or  any  general  excitation  that  may 
be  approximated  by  a piece-wise  linear  function. 


2.6.1  Steady  Excitation 

For  linear  systems  with  constant  properties  driven  by  a steady  excitation  the  response  of 
the  system  will,  eventually,  come  to  a steady  state  (i.e.,  {dT/dt}  = 0 ) given  by  the  solution  of: 

[K]{T}  = {E}  (2.18) 

2.6.2  Harmonic  Excitation 

For  linear  systems  with  constant  properties  driven  by  a steady  harmonic  excitation  of  the 
form: 

{E}  = Re(  {E*}  e'wt ) ; i = fT  (2.1 9a) 

where; 

{E*}  = the  complex  excitation  vector  (i.e.,  the  excitation  at  each  DOF  j,  E*j,  is 

represented  in  terms  of  amplitude  or  modulus,  E’j,  and  phase  angle  or  argument 
, <f>,  as; 

E*j  = E'j  e'4>  = E’j  (cos(<)>)  + i sin(<j>) ) (2.19b) 

to  = circular  frequency  of  excitation  [=]  radians/time 
the  response  of  the  system  will,  eventually,  come  to  a steady  harmonic  response: 

{T}  = Re(  (T* } e'wt ) (2.20) 
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given  by  the  solution  of: 


[K*]{T*J  = 

{E*} 

(2.21a) 

where: 

[K‘]  = [K  + itoC] 

(2.21b) 

m 

= complex  temperature  response 

= defined  in  terms  of  amplitude,  T'j , and  phase  angle  <|> , as: 

T*i 

= T'j  ei(l>  = Tj  (cos(4>)  + i sin(<t>) ) 

(2.21c) 

2.6.3  General  Excitation 

A finite  difference  scheme  for  the  approximate  integration  of  the  semi-discrete  equation 
(2.8)  may  be  developed  by  dividing  time  domain  into  discrete  steps: 

Wl  = tn  + 8t  ; n = 0,1 ,2,3  ...  (2.22) 

tQ  = initial  time 

where: 

5t  = integration  time  step  (often  constant  but  may  be  variable) 
demanding  the  satisfaction  of  equation  (2.8)  at  each  of  these  steps: 


[K]{T}n+i  + [C]{dT/dt}n+i  = (E}n+i 


(2.23) 


where: 

OWl  - (T(tn+i)} 

(dT/dt}n+i  = {dT(tn+i)/dt} 

(E}n+1  - (E(tn+1)} 


and  substituting  into  this  equation,  (2.23),  the  consistent  difference  approximation 
represented  by: 


fOn+1  - {T}n  + (1-a)5t{dT/dt}n  + a8t(dT/dt}n+1  (2.24) 

where: 

0 < a < 1 

a = 0 corresponds  to  the  Forward  Difference  scheme 
a = 1/2  corresponds  to  the  Crank-Nicholson  scheme 
a = 2/3  corresponds  to  the  Galerkin  scheme 
a = 1 corresponds  to  the  Backward  Difference  scheme 

a general  implicit  finite  difference  scheme  is  formulated: 


13 


[ a5t[K]  + [C]  ]{dT/dt}n+1  - {E}n+i  - [K]{  {T}n  + (1-o)5t{dT/dt}n  } 


(2.25a) 


or,  equivalently: 

[ [K]  + (1/a5t)[C]  ]{T}n+i  - {E}n+1  + (1/a5t)[C]{T}n  + ( 1 -cx)5t[C]{dT/dt}n  (2.25b) 

Computationally  it  is  strategic  to  implement  this  general  finite  difference  scheme,  equation 
(2.25),  as  a three  step  predictor-corrector  algorithm: 


rnn+1  = {Tjn  + (i-o)5t(dTdt}n 

; predictor  step 

(2.26a) 

[ a5t[K]  + [C]  ]{dT/dt}n+i  = {E}n+1  - [K]{T}n 

; (i.e.,  equation  (2.25a)  ) 

(2.26b) 

{T }n+i={T}n+i  +a5t  {dT/dt}n+1 

; corrector  step 

(2.26c) 

It  should  be  noted  that  this  algorithm  is  self-starting.  Given  initial  conditions,  (T(t0)}  , 
equation  (2.23)  may  be  solved  to  obtain  an  estimate  of  the  initial  rate  of  change  of  nodal 
temperatures,  {dT(t0)/dt}  , and  the  first  predictor  step,  equation  (2.26a),  may  then  be 
computed. 

This  predictor-corrector  scheme  has  been  analyzed  by  Taylor  [1]  and  Huebner  [5]  and  a 
more  general  predictor-multicorrector  scheme  that  includes  this  implicit  scheme  has  been 
analyzed  by  Hughes  [6].  For  a>  1/2  this  scheme  leads  to  an  unconditionally  stable 
(approximate)  solution;  for  a >3/4  (approximately)  leads  to  an  unconditionally  stable  non- 
oscillatory  solution;  beyond  this  Taylor  makes  some  recommendations  regarding  selection  of 
a and  step  size,  St,  to  limit  error  while  minimizing  computational  effort.  In  the  program 
DTAM1  the  default  value  of  a is  set  to  0.75,  and  may  be  reset  by  the  user,  and  an  estimate 
of  the  time  step  needed  to  limit  error  is  reported  (for  given  initial  conditions)  using  a method 
developed  by  Taylor  [1]. 


2.7  Specific  Element  Equations 

This  section  presents  the  element  equations  for: 

-ID  two-node  thermal  resistance  elements, 

- single-node  lumped  capacitance  elements, 

- two-node  fluid  flow  loop  elements, 

-ID  two-to-four  node  isoparametric  conduction  Finite  Elements, 

- 2D  four-node  isoparametric  conduction  Finite  Elements  (planar  and  axisymmetric),  and 

- variable-node  mean  radiant  temperature  element. 


2.7.1  Resistance  Elements 
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Resistance  elements  may  be  developed  directly  from  fundamental  heat  transfer  theory.  A 
2-node  resistance  element  results  that  is  defined  by  an  element  conductance  matrix  of  the 
general  form: 


[ k®]  = (A/R) 
where: 

A = the  area  available  for  heat  transfer 
R = the  element  resistance 


1 -1 
-1  1 


(2.27) 


Fia.  4 2-Node  Resistance  Element 

2.7.2  Lumped  Capacitance  Element 

Lumped  capacitance  elements  may  be  developed  directly  from  fundamental  heat  transfer 
theory.  A 1-node  capacitance  element  results  that  is  defined  by  an  element  capacitance 
matrix  of  the  general  form: 

[ce]  = (me)  [1]  (2.28) 

where: 

m = mass  of  element 

c = specific  heat  capacity 


j2  m , c 

Fia.  5 Lumped  Capacitance  Element 

2.7.3  Flow  Loop  Elements 

Energy  transferred  by  convective  flow  may  be  modeled  with  flow-related  elements.  Flow- 
related  element  equations  are  best  formulated  in  terms  of  nodal  temperatures  and  the  net 
rate  of  enthalpy  change,  henet_j  , (i.e.,  enthalpy  transported  from  the  node,  i,  into  the 
element).  For  the  simple  case  of  incompressible  fluid  flow  from  a single  node  to  another, 
ignoring  fluid  capacitance  effects,  the  element  equation  takes  the  form: 
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(2.29) 


where: 

w-|_2  = the  mass  flow  rate  from  node  1 to  node  2 

Cp  = the  specific  heat  capacity  at  constant  pressure  of  the  fluid 


h 


T 


e & 


net-1 


w, 


hnet-2 

t! 


O 


P'  '1-2 

Fia.  6 2-Node  Simple  Flow  Element 


In  principal,  element  equations  relating  rate  of  enthalpy  change  to  nodal  temperatures  may 
assembled,  along  with  nonflow  related  elements,  to  form  system  equations  governing 
coupled  flow  and  nonflow  heat  transfer.  Flow  related  element  equations  are,  however, 
characteristically  nonsymmetric,  as  in  the  case  presented  above.  The  program  DTAM1  is 
organized  to  take  advantage  of  the  property  of  symmetry  of  nonflow  elements  and, 
therefore,  even  this  simple  flow  related  element  may  not  be  added  to  its  library  of  elements. 
A subassemblage  of  two  simple  flow  elements,  creating  a flow  loop,  however,  results  in  a 
symmetric  system  of  equations  of  the  form: 


(2.30) 


where: 

w = the  mass  flow  rate  from  node  1 to  node  2 and  return,  (here  we  have  used 
superscripts  el  and  e2  indicate  the  two  simple  flow  elements  el  and  e2  of  the  subassembly) 


C , w 


Fig.  7 2-Node  Flow  Loop  Element 


The  flow  loop  element  equations  are  seen  to  be  similar  in  form  to  the  simple  resistance 
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element  equations  and,  therefore,  flow  loops  may  be  modeled  with  resistance  elements  by 
setting  the  quotient  (A/R)  equal  to  the  product  of  (wCp),  as  is  commonly  done  in  thermal 
network  analysis.  This  approach  to  formulating  flow  loop  equations,  based  upon  a 
consideration  of  rate  of  enthalpy  change,  leads  to  a entire  class  of  new  thermal  elements 
that  is  currently  under  investigation. 


2.7.4  Transient  Conduction  Finite  Elements 


The  solution  of  transient  heat  conduction  problems  using  the  finite  element  method  has 
been  discussed  by  several  authors  [1,  3,  5,  7,  8,  9,  & 10].  Over  a given  region,  Q , bounded 
by  a surface,  r , of  a solid  continuum,  the  spacial  and  temporal  variation  of  temperature, 
T(x,y,z,t)  , due  to  conduction  is  governed  by  the  so-called  transient  heat  conduction 
equation: 


(k  v— — 
ax  xdx 


dy 


(2.31a) 


with  temperature  boundary  conditions  specified  on  part  of  the  boundary,  r-|  : 


T = Tp  (x,  y,  z,  t=0)  ; on  r-j  , t > 0 


(2.31b) 


direct  heat  flow  and/or  convective  heat  flow  specified  on  part  of  the  boundary,  r2  : 

kxf^nx+  kVT~  nv+  kz  n z = q(x,y,z,t)  + h(x,y,z,t)  (T^-T)  ; onr2,t  > 0 

dx  ydy  y *dz  (2.31c) 


and  initial  conditions: 


T = T0(x,  y,  z,  t=0)  ; in  Q 
where; 


(2.31d) 


kx  » ky  • kz 

P 

C 


q* 

TP 

nx  , ny  , nz 

q 

h 

Too 


= the  principal  conductivities  for  a thermally  anisotropic  continuum 
= the  density 

= the  specific  heat  capacity 
= the  internal  heat  generation  rate  per  unit  volume 
= the  prescribed  temperature  boundary  condition 
= the  direction  cosines  of  a normal  to  the  surface  at  (x,  y,  z) 

= the  prescribed  heat  flow  rate  boundary  condition  on  r2 
= the  convective  heat  transfer  coefficient  on  r2 
= the  temperature  of  the  external  environment  that  drives  the 
convective  heat  flow 


In  the  finite  element  method  a solution  for  the  temperature  field,  T(x,  y,  z,  t),  is 
approximated  by  assumed  functions  that  are  defined  over  a small,  but  finite,  part  of  the 
region  (i.e.,  the  finite  element).  These  functions  have  the  general  form: 
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T(x,  y,  z,  t)  = {N(x,  y,  z)}{Te(t)} 


(2.32) 


where  the  spacial  {unctions,  {N(x,  y,  z)},  - the  so-call  shape  functions  - are  selected  to 
satisfy  continuity  requirements  across  inter-element  boundaries.  The  spacial  variation  of  the 
temperature  field  is  thus  approximated  in  a piecewise  manner,  using  the  shape  functions,  and 

the  temporal  variation  is  captured  by  the  variation  of  discrete  values  of  temperature,  {T®}, 
that  correspond  to  the  element  temperature  DOFs  discussed  above. 

Using  the  approximation  represented  by  equation  (2.32)  and  the  Galerkin  variation  of  the 
method  of  weighted  residuals,  the  transient  heat  conduction  equation  is  transformed  to 
equations  of  the  form  of  (2.8)  that  may  be  assembled,  in  the  manner  discussed  above,  from 
element  equations  of  the  form  of  equation  (2.5)  where,  now: 

- the  element  conductance  matrix  coefficients  for  elements  in  Q.  and  on  T-j  are: 


(2.33a) 


- the  element  conductance  matrix  coefficients  for  elements  on  r2  are: 


(2.33b) 


- the  element  capacitance  matrix  coefficients  for  all  elements  are: 


(2.33c) 


the  element  internal  generation  rate  vector  components  are: 


(2.33d) 


- the  system  direct  heat  flow  rate  vector  is  augmented  by  the  components: 


(2.33e) 
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where  Qj®  , the  element  contribution  to  the  system  direct  heat  flow  rate  at  element  node  i, 
includes  two  terms  - a direct  surface  gain  term  and  a convective  gain  term,  respectively. 

A practically  unlimited  variety  of  specific  conduction  elements  may  be  formulated  using 
these  equations;  no  restrictions  on  the  number  of  nodes  per  element  have  been  made  and 
the  compatibility  requirements  placed  upon  admissible  shape  functions  is  not  very  restrictive. 
Furthermore,  even  after  decisions  regarding  the  number  of  nodes  per  element  and  the 
choice  of  shape  functions  has  been  made  the  analyst  may  choose  to  employ  analytical  or 
various  numerical  integration  strategies  to  evaluate  element  matrices  thereby  generating 
subclasses  of  element  types. 

Elements  having  distorted  geometries  and/or  having  nonuniform  node  spacings  may  be 
conveniently  formulated  by  mapping  the  global  (i.e.,  actual  or  physical)  coordinate  geometry 
(x,y)  to  an  undistorted  local  geometry  (r,s)  and  evaluating  the  integral  expressions,  above, 
in  the  undistorted  geometry  space  using  relatively  simple  shape  functions  defined  in  the 
undistorted  geometry  space.  If  the  mapping  from  the  global  geometry  to  the  local  geometry 
is  defined  using  the  same  shape  functions  used  to  approximate  the  temperature  field  (i.e.,  in 
equation  (2.32))  then  the  element  formulation  is  said  to  be  isoparametric.  (An  introduction 
to  isoparametric  element  formulation  may  be  found  in  [8]  pp. 193-230  and  the  practical 
implementation  of  the  approach  is  presented  in  [11].) 

Two  finite  elements  are  provided  in  the  program  DTAM1 ; a variable  2 to  4-node 
isoparametric  one-dimensional  conduction  finite  element  and  a 4-node  isoparametric  two- 
dimensional  finite  element. 

Two  to  Four-Node  One-dimensional  Isoparametric  Elements 

A variable  two  to  four-node  isoparametric  conduction  finite  element  may  be  formulated  for 
one-dimensional  heat  transfer  by  using  the  following  shape  functions,  {N}  = {N-|,  N2,  N3,  N4}, 
defined  relative  to  a biunit  element  in  local  coordinate  r as  (see  [3]  page  199); 


Include  only  if  Include  only  if  node 

node  3 is  present  3 and  4 are  present 


N,  = ~r(1  - r) 

■ “0  - r2) 

2 

2 

N2  = ^(1  - r) 

-l(1-r2) 

N3  = (1  - r2) 

+ -^-(-9r3  + r2  + 9r -1)  (2.34) 

1 6 

* — (9r3  + r2  - 9r  - 1) 

1 6 

♦ (27r3  + 7r2  - 27r  - 7) 

1 6 


N4  = — ( - 2 7 r 3 - 9r2  + 27r  + 9) 
1 6 


where; 

r = the  local  coordinate  defined  as  indicated  below: 
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Global  Coordinates 


Local  Coordinates 


3 


X) 


r = - 1 + 1 

r = - 1 0 + 1 

r = - 1 -v3  +1/3  + 1 


2- node 

3- node 

4- node 


Fig.  8 Two  to  Four-Node  Element  Coordinate  Systems 

In  the  isoparametric  formulation  the  element  global  coordinate,  x,  is  related  to  the  element 
local  coordinate,  r,  through  a transformation  defined  by  the  shape  functions  given  above: 

x = {N}{x}  (2.35) 

where: 

{N}  = {Ni  , N2}  for  two-node  elements 

= {Ni  , N2  , N3}  for  three-node  elements 
= {Ni  , N2  , N3  , N4}  for  four-node  elements 
{x}  = {xi  , X2}  for  two-node  elements 

= {xi  , X2  , X3}  for  three-node  elements 
= {xi  , X2  , X3  , X4}  for  four-node  elements 

The  integral  equations  (2.33)  may  then  be  rewritten  in  terms  of  the  local  coordinates  as: 


rx2 

dN,  dNj 
I dx  dx 


Cy  = A pc 


/X2  r +1 

NjNjdx  = Ape  J NjNj 

A 1 "'-I 


J dr 


(2.36a) 


(2.36b) 


(2.36c) 


where: 

A = the  area  available  for  heat  transfer 
k = kx  ; the  element  conductivity 

J = dx/dr  = d/dr{N}{x}  ; the  Jacobian  of  the  coordinate  transformation. 
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It  is  convenient  to  evaluate  these  integrals  numerically.  Gauss  quadrature  may  be  used  to 
exactly  evaluate  these  integrals,  but  the  resulting  equations  will  tend  to  be  "stiff"  and,  as  a 
result,  will  tend  to  demonstrate  artificial  oscillatory  behavior  [12,13].  Hughes  has  shown 
that  this  problem  can  be  mitigated  through  approximate  evaluation  of  the  capacitance  matrix 
using  appropriate  quadrature  rules  [6].  A diagonal  or  lumped  capacitance  matrix,  which  tends 
to  inhibit  artificial  oscillatory  behavior  at  the  expense  of  accuracy,  will  result  if  quadrature 
points  for  the  evaluation  of  the  capacitance  matrix,  (2.36b),  are  selected  to  be  equal  to  the 
element  nodal  coordinates  in  the  local  coordinate  system. 

For  the  two-node  element  the  element  matrices  include: 

- the  element  conductance  transfer  matrix: 

[k6)  = (Ak/L) 


1 -1 
-1  1 


12.37) 


- the  element  thermal  capacitance  matrix: 


[c®|  = (ApcL) 


(1/2 


r)  r 
(1/2 - r) 


(2.38) 


where: 

r = 0 corresponds  to  a lumped  capacitance  model  that  results  from  using  the 
quadrature  points  (-1  , +1) 

= 1/12  corresponds  to  a higher  order  capacitance  model  that  results  from 
using  the  quadrature  points  (-V  2/3  , +Y2/3  ) 

= 1/6  corresponds  to  a consistent  capacitance  model  that  results  from  using 
Gauss  quadrature  (i.e.,  quadrature  points  -1/f3  , +1/VT  ) to  exactly  evaluate 
equation  (2.36b) 
and 


- the  element  internally  generated  heat  flow  rate  vector  (e.g.,  a radiant  floor)  is: 

(q8)  = q‘L/2^  j 


(2.39) 


It  may  be  shown  that  the  higher  order  capacitance  matrix  will  minimize  error  in  the 
evaluation  of  temperature  [6]. 


Ti  A , L , k j, 
qf  P > c , q*  qf 

Fia.  9 One-Dimensional  Two-Node  Conduction  Finite  Element 
When  used  with  the  implicit  numerical  integration  scheme  provided  in  the  program 
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DTAM1,  the  higher  order  capacitance  model,  r = 1/12,  results  in  the  greatest  accuracy  for  a 
given  time  step.  The  choice  of  r=0  results  in  equations  identical  to  those  associated  with  the 
finite  difference  method  solutions. 


Four-Node  Two-Dimensional  Isoparametric  Element 

Element  matrices  for  a four-node  two-dimensional  isoparametric  element  may  be  evaluated 
in  a similar  manner  as  that  presented  above  for  the  one-dimensional  isoparametric  element. 
The  program  DTAM1  provides  planar  and  axisymmetric  elements  that  are  evaluated  using 
Gauss  integration  of  integral  expressions  corresponding  to  equations  (2.33).  As  Gauss 
integration  results  in  the  exact  evaluation  of  these  integrals  these  elements  will  result  in 
"stiff"  equations,  which  tend  to  demonstrate  artificial  oscillatory  behavior  in  regions  where 
the  temperature  field  has  rapidly  changing  and  pronounced  gradients. 


2.7.5  Variable  Node  Mean  Radiant  Temperature  Elements 

A variable  node  nonlinear  radiative  element  may  be  derived  directly  from  Joseph  Carroll's 
mean  radiant  temperature  network  method  [14,  15,  16]  that  has  a single  node  for  each 
surface  considered  and  an  additional  node: 

"The  fictitious  Mean  Radiant  Temperature  node  in  the  center  of  the  network  [that]  acts 
as  a clearinghouse  for  all  radiative  interchanges." 

The  surface  nodes  are  linked  to  the  MRT  node  by  simple  resistance  elements  whose 
resistances  are  defined,  nonlinearly,  by  Carroll's  method.  The  element  conductance  transfer 
matrix  for  this  variable  node  radiative  element  will  have  the  form: 


V 


Planar  Element 


Axisymmetric  Element 


Fig.  10  Four-Node  Isoparametric  Finite  Element 
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[0  = hb 


KAjF’j)  (-A  -|  F ’ 1 ) (-A2  F '2 ) 
(“A i F ' ) (A-jF’t)  0 

(~A 2 F '2 ) 0 (A2  F '2 ) 


(-AnF'n)  0 


0 


("AnF’n) 

0 

0 


(An  F n ) 


(2.40) 


where: 

The  element  DOF  are  numbered  i=0, 1,2,3  ...  n with  DOF  i=0  corresponding  to  the  MRT  node. 


hb  = 4a  T3 

a = Stefan-Boltzmann  Constant 
T = the  average  absolute  temperature  of  surface 
Aj  = area  available  for  radiative  exchange  from  surface  i 

F'j  = radiant  interchange  factor  between  surface  i and  the  MRT  node 

= 1/ ( 1/Fj  + (1*Ej)  / £j) 


Fj  = 1/(  1—  A j Fj/  X(AjFj)  ) ; determined  iteratively 
£j  = emissivity  of  surface  i 


Fia.  11  Variable  Node  MRT  Element 
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3.  Summary  and  Conclusions 


An  approach  to  building  energy  simulation  has  been  presented  that  is  based  upon  discrete 
analysis  techniques  used  in  other  fields  of  physical  system  simulation  wherein  system 
behavior  is  governed  by  an  overriding  equilibrium  principle  (e.g.,  mass,  momentum,  or  energy 
conservation)  that  is  used  to  define  an  assembly  process  that  allows  the  formation  of 
system  equations  from  individual  element  equations.  The  element  equations  are  developed  to 
model  the  behavior  of  discrete  regions  (e.g.,  finite  elements)  or  discrete  processes  within 
the  system.  In  the  present  context  the  governing  equilibrium  principle  is  the  conservation  of 
energy  and  element  equations  have  been  presented  that  model  conductive,  convective,  and 
radiative  energy  transport  and  that  account  for  thermal  storage  processes  in  building 
systems. 

The  discrete  thermal  analysis  theory  and  its  practical  implementation  in  the  program 
DTAM1  has  been  presented  to  provide  a first  demonstration  of  this  very  general  approach. 
It  has  been  shown  that  element  equations  may  be  developed  from  very  basic  considerations 
(e.g.,  the  ID  thermal  resistance  element  and  the  simple  capacitance  element),  from  Finite 
Element  approximations  to  governing  partial  differential  equations  (e.g.,  the  ID  and  2D 
isoparametric  conduction  Finite  Elements),  and  from  more  heuristic,  ad  hoc  physical 
arguments  (i.e.,  the  variable-node  MRT  element).  Although  the  discussion  was  limited  to 
element  formulations  with  constant  thermal  properties  the  extension  to  nonconstant  and 
nonlinear  thermal  properties  is  straightforward.  Consequently,  a practically  unlimited  variety 
of  elements  could,  conceivably  be  developed  following  the  examples  presented  in  this  report. 
In  particular  the  nonsymmetric  flow  element  discussed  in  section  2.7.3  should  lead,  quite 
naturally,  to  the  development  of  a variety  of  HVAC  component  elements. 

One  may  also  approach  the  formation  of  system  equations  based  upon  either  conduction 
response  function  theory  or  transmission  matrix  theory  from  an  element  assembly  approach 
and,  thus,  consider  the  development  of  response  function  elements  and/or  transmission 
matrix  elements.  In  these  cases,  however,  the  system  equations  that  result  will  be  algebraic, 
rather  than  differential,  the  extension  to  nonlinear  situations  is  not  so  straightforward,  and 
mixed  assemblies  of  elements  based  upon  the  approach  presented  here  with  those  based 
upon  these  heat  transfer  theories  will  only  be  possible  in  special  cases.  Nevertheless,  by 
taking  this  approach  to  the  formulation  of  these  approaches  one  is  able  to  bring  the  full 
range  of  building  energy  simulation  possibilities  under  one  unifying  theory. 

Steady  state,  steady  periodic,  and  transient  analysis  options  for  thermal  systems  with 
constant  thermal  properties  (i.e.,  linear  thermal  systems)  have  been  discussed.  Additional 
solution  options  may  also  be  considered  for  linear  thermal  systems  and  a variety  of 
techniques  exists  for  the  analysis  of  thermal  system  with  nonconstant  and/or  nonlinear 
thermal  properties;  Table  1 outlines  some  of  these  possibilities. 

From  a theoretical  point  of  view  the  discrete  thermal  analysis  method  serves  to  unify 
existing  building  energy  simulation  theory  and  place  it  on  a rigorous  base.  In  so  doing  the 
building  energy  simulation  community  will  be  in  a better  position  to: 

- compare  advantages  and  disadvantages  of  existing  building  energy  simulation 
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methods, 

- make  use  of  the  enormous  volume  of  research  literature  presently  available  in  the 
general  area  of  discrete  analysis  techniques  (e.g.,  the  Finite  Element  and  related 
numerical  methods),  and 

- move  on  to  the  more  demanding  task  of  simulating  coupled  thermal,  flow,  and  HVAC 
system  behavior  for  indoor  air  quality  analysis. 

From  a computational  point  of  view  the  discrete  thermal  analysis  method,  being  based 
upon  distinct  formal  operations  (e.g.,  formation  of  element  matrices,  assembly  of  element 
matrices,  solution  of  system  equations,  & post  evaluation  of  element  results)  using  a library 
of  thermal  elements,  leads  to  highly  modular  program  structures  and  could  be  the  basis  of  a 
higher  level  programming  language  for  building  energy  simulation  program  development. 

From  a practical  point  of  view  the  discrete  thermal  analysis  method  allows  an  analyst  to 
make  use  of  all  the  "tools"  that  have  emerged  from  building  energy  simulation  research 
community,  to  create  detailed  or  simple  idealization  - as  suits  the  purpose  - that  may  readily 
be  modified  as  design  evolves  (i.e.,  through  the  addition,  modification,  or  deletion  of 
elements),  and,  for  a given  idealization,  to  consider  simple  to  complex  models  of  excitation 
idealization  as  appropriate  to  the  stage  of  design  or  level  of  analytical  inquiry. 


Excitation  Model 


Table  1 Solution  Options  for  The  Governing  Discrete  Thermal  Analysis  System  Equations 

t K ] {T } + [ C]{dT/dt } = {E} 
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DTAM1  USERS  MANUAL 


I.  General  Instructions 


Thermal  analysis  of  a building  system,  using  DTAM1,  involves  three  basic  steps: 


Step  1 Idealization  of  the  Building  System  and  Excitation 


Fig.  1.1  Idealization  of  the  Building  System  and  Excitation 


Idealization  of  the  building  thermal  system  involves  definition  of  a coordinate  system, 
discretization  of  the  system  as  an  assemblage  of  appropriate  thermal  elements  connected  at 
system  nodes,  identification  of  boundary  conditions,  numbering  of  system  nodes  optimally 
(i.e.,  to  minimize  the  bandwidth  of  system  equations),  and  estimation  of  the  bandwidth  of 
the  system  equations  (i.e.,  the  largest  coupled  DOF  number  difference  plus  1).  Selection  of 
appropriate  thermal  elements  involves  some  skill  and  considerable  knowledge  of  the 
characteristics  of  building  thermal  behavior  and  the  numerical  and  behavioral  characteristics 
of  specific  thermal  elements. 

The  excitation  - specified  temperatures,  direct  nodal  heat  flow  rates,  and  external 
environmental  temperatures  associated  with  convection  - may  be  modeled  to  be  steady,  to 
vary  harmonically,  or  defined  in  terms  of  arbitrary  piecewise-linear  time  histories.  For  the 
latter  case  initial  conditions  of  nodal  temperatures  will  have  to  also  be  specified.  (For  those 
problems  where  these  initial  conditions  are  not  known  apriori  the  analyst  may  use  a false 
start-up  period  to  attempt  to  determine  initial  conditions.) 

Step  2:  Preparation  of  Input  Command/Data  File 


The  program  DTAM1  is,  in  essence,  a command  processor.  The  program  reads  an  ASCII 
text  file  that  contains  commands  and  associated  data,  collected  together  in  distinct  data 
groups,  that  define  the  thermal  idealization.  The  command/data  input  file  may  be  prepared 
with  any  available  ASCII  text  editing  program  and  given  a file  name,  <filename>,  specified  by 
the  user.  The  <filename>  must,  however,  consist  of  8 or  less  alphanumeric  characters  and 


Flq.  1.2  Preparation  of  Input  Command/Data  File 
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can  not  include  an  extension  (i.e.,  characters  separated  from  the  filename  by  a period, ). 


Step  3 : Execution  of  DTAM1 


ooioft 
01100 
00011  I 

<filename>.DBS 


ooio^i 

01100 

00011 


<filename  >.TMP 


<filename> 


<filename>.OUT 


Fio.  1.3  Execution  of  DTAM1 

DTAM1  is  then  executed.  It  will  request  the  name  of  the  command/data  input  file  and 
proceed  to  form  element  and  system  arrays  and  compute  the  solution  to  the  posed  problem. 
DTAM1  reads  the  ASCII  input  file  and  creates  a single  ASCII  (i.e.,  printable)  output  file  and 
up  to  four  binary  files: 

Files  Read 

<filename>  an  ASCII  input  file  specified  by  the  user  that  contains  problem  data, 

Files  Created 

<filename>.OUT  a printable  ASCII  output  file  that  contains  analysis  results, 

<filename>.DBS  a binary  file  containing  system  arrays  that  may  be  read  by  CAL-80  program 
segments  for  further  processing  of  system  arrays, 

<filename>.TMP  a binary  file  containing  the  computed  time  histories  of  nodal  temperatures, 
<filename>.PRE  a binary  file  used  for  temporary  disk  storage  of  system  arrays  used  by  the 
predictor-corrector  integration  subroutines  in  DTAM1 , and 
<filename>.ELM  a binary  file  used  for  temporary  disk  storage  of  element  arrays. 

The  results  of  an  analysis,  <filename>.OUT,  may  be  conveniently  reviewed  using  an  ASCII 
editor  and,  from  the  editor,  portions  or  all  of  the  results  may  be  printed  out. 
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II.  Preparation  of  Input  File 


Input  data  is  defined  in  terms  of  data  groups  that  are  delimited  by  a data  group  separator 
label,  that  may  be  thought  to  be  a command,  and  an  end  flag.  With  the  exception  of  the 
problem  initiation  data  group,  data  groups  may  be  ordered  as  the  user  sees  fit. 
Command/data  input  files  follow  the  free-field  input  conventions  of  the  CAL-SAP  software 
development  system  (see  Appendix  A). 

The  specific  data  groups  needed  to  be  included  in  a command/data  input  file  will  depend 
upon  the  details  of  the  thermal  idealization,  the  nature  of  the  excitation,  and  the  type  of 
analysis  to  be  computed.  For  the  dynamic  analysis,  using  the  predictor-corrector  integration 
option,  an  input  file  will  have  the  following  general  form  (where  nn  is  a data  value): 


Column 

Problem  Initiation: 


1 

DTAM1 

<problem  labeling  information> 

NODEmn  B AND  = nn 
SAVE 


Constraint  Data  Group:  CONSTRAIN 

nn,nn,nn  M=nn 

END 


Coordinate  Data  Group:  COORDINATES 

nn  X=nn  Y=nn  Z=nn 

nn  X=nn  Y=nn  Z=nn  GEN=nn,nn,nn 

END 

Element  Data  Groups: 

Resistance  Elements:  RESIST 

nn  l=nn,nn  R=nn  A=nn 
nn  l=nn,nn  R=nn  A=nn  GEN=nn 

END 


Capacitance  Elements: 


Isoparametric  ID 
2-Node  Elements: 


Isoparametric  2D 
4-Node  Elements: 


CAPACI 

nn  l=nn  M=nn  C=nn 

nn  l=nn  M=nn  C=nn  GEN=nn 

END 

ISOI  D 

nn  l=nn,nn  K=nn  P=nn  C=nn  A=nn 
nn  l=nn,nn  K=nn  P=nn  C=nn  A=nn  GEN=nn 

END 

IS02D4 

nn  l=nn,nn,nn,nn  K=nn,nn,nn  C=nn  P=nn 
nn  l=nn,nn,nn,nn  K=nn,nn,nn  C=nn  P=nn 

END 


Boundary  Conditions  Group:  BOUNDARY  COND 
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nn.nn.nn  BC=nn 
nn,nn,nn  BC=nn 

END 

Convection  Boundary 
Conditions  Group: 

CONVECTION 
nn  l=nn,nn  H=nn 
nn  l=nn,nn  H=nn  GEN=nn 

Solution  Control 
- Predict-Correct.  Integ. 
Solution  Control  Data: 

END 

PREDICT 

TIME=nn,nn,nn  ALPHA=nn  NEFN  = nn 
nn.nn.nn  T=nn  Q=nn  TEX=nn 
nn.nn.nn  T=nn  Q=nn  TEX=nn 

END 

Initial  Conditions  Group: 

INITIAL  COND 
nn.nn.nn  T=nn 
nn.nn.nn  T=nn 

END 

Excitation  Function 
Data: 

EXCITATION 

T=nn  EFN=nn,nn,nn,  ... 

T=nn  EFN=nn,nn,nn,  ... 

Report  Control  Data: 

END 

REPORT 
WRITE  PINT=nn 
nn.nn.nn 

END 

Details  are  given  on  the  following  pages  for  each  data  group. 
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A.  Problem  Initiation 


An  analysis  problem  must  be  initiated  by  the  following  lines  of  data: 

DTAM1 

<problem  labeling  information>  : 

N0DE=I1  BAND  = 82 

[SAVE] 

where: 

11  = the  number  of  nodes 

12  = the  (half)  bandwidth;  the  largest  DOF  number  difference  between  coupled 
nodes  plus  1 

If  the  keyword  SAVE  is  found  at  the  beginning  of  the  fourth  line,  system  arrays  will  be 
written  to  the  binary  file  <filename>.DBS  that  may  be  read  by  CAL-80  program  segments 
for  further  processing.  CAL-80  |16]  is  a high-level  language  developed  to  be  used  as  an 
educational/research  tool  that  facilitates  basic  matrix  algebraic  operations  and  numerical 
analysis  methods  (e.g.,  equation  solving,  eigenanalysis,  fast  fourier  transformation,  etc.). 
CAL-80  may  be  obtained  from  the  National  Information  Service  for  Earthquake  Engineering, 
NISEE/Computer  Applications,  Davis  Hall,  University  of  California,  Berkeley,  California 
94720,  (415)  642-5113. 


b.  Const raJinl-Jaia  - Group 

System  temperature  DOFs  may  be  constrained  to  be  equal  (e.g.,  see  Fig.  1)  by  the 
following  lines  of  data: 

CONSTRAIN 
11,12,13  M=S4 

END 

where: 

11,  12,  13  = first  node,  last  node,  node  number  increment  of  a group  of  slave 

nodes  constrained  to  be  equal  to  the  same  master  node;  13=1  default 
14  = the  master  node  to  which  the  other,  slave  , nodes  are  constrained 

Constraining  system  temperature  DOF  provides  one  means  of  using  one-dimensional 
elements  in  conjunction  with  two-dimensional  elements.  If  not  done  with  care,  constraining 
may  have  the  effect  of  destroying  the  compactness  of  system  arrays  and  thus  result  in 
unnecessary  computational  effort  in  solving  the  system  equations. 
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C.  Coordinate  Data  Group 


Coordinate  data  is  required  for  some,  but  not  all,  elements.  It  is  defined  by  the  following 
lines  of  data: 

COORDINATES 

II  X=R1  Y=R2  Z=R3  GEN=I2,I3,I4 
END 

where; 

II  = the  node  number 

R1  = the  x coordinate 

R2  = the  y coordinate 

R3  = the  z coordinate 

12,13,14  = first  node,  last  node,  node  number  increment;  14=1  default 

If  the  generation  data  is  given,  nodal  coordinates  are  generated  at  equal  intervals  along  a 
straight  line  between  the  first  and  last  nodes,  12  & 13,  incrementing  node  numbers  by  the 
generation  increment,  14.  Nodal  coordinates  not  given  will  be  assumed  to  be  zero. 
Coordinate  data  is  not  needed  for  all  elements;  see  specific  element  input  description  below 
for  details. 


Either  rectangular  and  cylindrical  coordinate  data  may  be  defined  but  right  hand 
coordinate  systems  must  be  used: 


X X (R) 


D.  Element  Data  Groups 

Element  data  is  entered  by  element  type.  Details  for  each  element  type  is  given  below. 
Each  element  type  follows,  however,  the  general  form: 

<element  name> 

II  1=12,13,  ...  GEN=I4  <element  properties>  [0=C1] 

END 
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where: 

II  = the  element  number 

12,13,...  = the  element  node  numbers 

14  = generation  increment  (default  = 1) 

Cl  = Y if  element  nodal  heat  flow  rate  results  are  to  reported  (i.e.,  output) 

= N if  element  nodal  heat  flow  rate  results  are  not  to  be  reported;  N default. 

Element  data  must  be  supplied  in  numerical  order.  Omitted  element  data  is  automatically 
generated  by  incrementing  the  preceding  node  numbers  by  the  current  generation 
increment.  Properties  of  generated  elements  will  be  set  equal  to  that  of  the  current  element. 

***  The  0=Y  option  is  presently  not  implemented.  *** 


1)  Resistance  Elements 

Conventional  2-node  resistance  element  equations  may  be  added  to  the  system  equations 
with  the  following  lines  of  data: 


RESIST 

II  1=12,13  GEN=I4  R=R1  A=R2 


END 


where: 

II 

12,  13 
14 
R1 
R2 


- the  element  number 
= the  element  node  numbers 
= generation  increment  (default  = 1) 
= the  element  resistance 
= the  area  available  for  heat  transfer 


2}  Lumped  Capacitance  Elements 

Conventional  1-node  capacitance  element  equations  (i.e.,  lumped  capacitance)  may  be 
added  to  the  system  equations  with  the  following  lines  of  data: 

CAPAC! 

SI  3=12  GEN=I3  M=R1  C=R2 
END 

where; 

S 1 = the  element  number 

12  = the  element  node  number 

13  = generation  increment  (default  = 1) 
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R1  = the  thermal  mass 

R2  = the  specific  heat  capacity  of  the  thermal  mass 


3)  Fluid  Flow  Loop  Elements 

Elements  equations  that  model  heat  transfer  via  a fluid  flow  loop  between  any  two  nodes 
may  be  added  to  the  system  equations  with  the  following  lines  of  data: 

FLOWLOOP 

II  1=12,13  GEN=I4  W=R1  C=R2 
END 

where: 

II  = the  element  number 

12,  13  = the  element  node  numbers 

14  = generation  increment  (default  = 1) 

R1  = the  mass  flow  rate  from  first  node  to  second  & second  node  to  first 
R2  = the  fluid  specific  heat  capacity  (at  constant  pressure) 


4)  Isoparametric  ID  Two  to  Four-node  Conduction  Finite  Elements 

Isoparametric  ID  two  to  four-node  conduction  finite  element  equations  may  be  added  to 
the  system  equations  with  the  following  lines  of  data: 

ISOID 

II  1=12, 13, [14, 15]  GEN=I6  X=R1  ,R2,[R3,R4]  K=R5  P=R6  C=R7  A=R8  [Q=R9] 

[R=— ] 


END 


where: 


II  = the  element  number 

12,13  = the  element  end  node  numbers  (for  node  number  order  see  sketch  below) 

14,15  «=  the  intermediate  nodes  (for  node  number  order  see  sketch  below): 

if  13  < 0 a two-node  linear  element  will  be  generated 

if  13  > 0 & 14  < 0 a three-node  quadratic  element  will  be  generated 

if  13  > 0 & 14  > 0 a four-node  cubic  element  will  be  generated 

!6  = generation  increment  (default  = 1) 

R1,R2  = the  element  end-node  coordinates  (i.e.,  along  the  length  of  the  element) 

R3,R4  = the  element  intermediate  node  coordinates 

R5  = the  element  conductivity 

R6  = the  element  weight  density 

R7  = the  element  specific  heat  capacity 

R8  = the  area  available  for  heat  transfer 
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R9  = the  internal  heat  generation  rate  per  unit  length  of  element 
R=...  = GAUSS  or  = RIO,  R11,  [R12.R13] 

- if  GAUSS  Gauss  quadrature  will  be  used  to  evaluate  the  element 
capacitance  matrix 

- RIO,  R11,  [R12.R13]  are  quadrature  points  to  be  used  to  evaluate 
element  capacitance  matrix;  the  number  of  quadrature  points  used  will  be 
equal  to  the  number  of  element  nodes  (this  assures  exact  integration  when 
Gauss  quadrature  points  are  used);  the  following  quadrature  points  will 
result  in  a lumped  capacitance  modeling: 


Element 
2-Node  Linear 

-1.00 

RIO 

1.00 

Ell 

R 1 2 

R12 

3-Node  Quadratic 

-1.00 

0.00 

1.00 

default 

4-Node  Cubic 

-1.00 

-1/3 

+1/3 

1.00 

default 

- for  2-node  linear  elements  the  quadrature  points  R = - V 2/3,  + Y2/3  will 
minimize  error;  these  quadrature  points  are  the  default  quadrature  points 
for  the  2-node  linear  element. 

1 2 

% % Two-Node  Element 

1 3 2 

^ Three-Node  Element 

1 3 4 2 

% % % % Four-Node  Element 

Fia.  11.2  One-Dimensional  Isoparametric  Element  Node  Numbering 

NOTE: 

a)  One  may  implement,  in  effect,  a ID  finite  difference  model,  corresponding  to  a 
linear  difference  approximation  over  the  interval  of  the  element  length,  by  using 
lumped  capacitance  two-node  elements. 

b)  When  using  3-node  quadratic  and  4-node  cubic  elements  lumped  capacitance 
modeling  may  provide  more  reliable  results  than  Gauss  quadrature  capacitance 
modeling  (see  Example  B1  discussed  in  Appendix  B). 


5)  Isoparametric  2D  Four-node  Conduction  Finite  Elements 

Isoparametric  2D  4-node  conduction  finite  element  equations  may  be  added  to  the  system 
equations  with  the  following  lines  of  data: 

IS02D4 

II  1=12,13,14,15  GEN=I6  K=R1,R2,R3  C=R4  P=R5  [Q=R6]  [D=R7]  [REG=C1] 
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END 


where; 

II  = the  element  number 

12,13,14,15  = the  element  node  numbers  ordered  counterclockwise  around  the 

element 

R1.R2.R3  = the  components  of  the  conductivity  tensor,  Kxx,Kyy,KXy;  for 
homogeneous,  isotropic  materials  Kxx=Kyy=k  , Kxy=0;  for  anisotropic  materials 
aligned  to  the  chosen  coordinate  system  kxx  = kx  , kyy  = ky  , kxy  = 0 

R4  = the  specific  heat  capacity 

R5  = the  weight  density 

R6  = the  internal  heat  generation  per  unit  volume:  R6=0.0  default 

R7  = the  thickness  (planar  regions)  or  subtended  angle  (axisymmetric  regions)  of 
the  element:  R7=1 .0  default 

Cl  = PLN  for  planar  regions  or  AXI  for  axisymmetric  regions:  PLN  default. 

Triangular  elements  are  possible:  they  are  identified  by  repeating  any  two  nodal  point 
numbers  (e.g.,  I,J,K,K). 


E.  Boundary  Conditions  Data  Group 

Boundary  condition  flags  are  established  with  the  following  lines  of  data: 

BOUNDARY  CONDITIONS 
11,12,13  BC=C1 

END 

where: 

11,12,13  = the  first  node,  last  node,  node  number  increment  of  a series  of  nodes  with 
identical  boundary  conditions 

Cl  = Q for  heat  flow  rate  prescribed  boundaries:  T for  temperature  prescribed 
boundaries;  Q default 

The  heat  flow  rate  at  the  temperature  - but  not  both  - may  be  specified  at  each  node  to 
establish  boundary  conditions  of  prescribed  temperature  or  heat  flow  rate.  These  conditions 
are  specified  with  as  many  lines  as  needed. 

If  this  data  group  is  omitted  all  nodes  will  be  assumed  to  be  heat  flow  rate  DOF. 


F.  Convection  Boundary  Conditions  Data  Group 
Convection  boundary  conditions  are  established  by  the  following  lines  of  data: 
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CONVECTION  BOUNDARY  CONDITIONS 
II  1=12,13  GEN=I4  H=R1  [D=R2]  [REG=C1] 

END 


where: 

II 

12,13 

14 

R1 

R2 

Cl 


= the  convective  surface  segment  number 
= the  surface  segment  end  nodes 
= the  generation  increment:  14=1  default 
= the  convective  transfer  coefficient 

= the  the  thickness  (planar  regions)  or  subtended  angle  (axisymmetric  regions) 
of  the  element:  R2=1.0  default 

= PLN  for  planar  regions  or  AXI  for  axisymmetric  regions:  PLN  default. 


Convective  boundary  conditions  may  be  specified  for  surface  segments  defined  by  pairs 
of  nodes  that  have  been  specified  to  be  part  of  the  Q boundary.  For  these  nodes,  then,  it 
is  possible  to  account  for  both  general  direct  heat  flow  excitation  and  convective  excitation 
(e.g.,  direct  solar-to-surface  gain  and  convective  air-to-surface  gain).  These  conditions  are 
specified  with  as  many  lines  as  needed. 


(Convective  boundary  conditions  would  reasonably  only  be  associated  with 
subassemblages  of  two-dimesional  conduction  elements.  Furthermore,  one  would  reasonably 
only  specify  plane  convective  transport  for  subassemblages  of  plane  elements  and 
axisymmetric  convective  transport  for  subassemblages  of  axisymmetric  elements.  No 
attempt  is  made,  however,  to  check  the  consistency  of  the  model.) 

G.  Solution  Control 


Data  groups  needed  to  control  the  solution  procedure  and  define  the  excitation  are 
explained  below. 

1)  Steady  State  Excitation/Response 

The  response  of  the  system  to  steady  excitation  may  be  computed  by  including  the 
following  lines  of  data: 


STEADY 

11,12,13  T=R1  Q=R2  TEX=R3 


END 


where: 

11,12,13  * the  first  node  number,  last  node  number,  and  node  number  Increment  of  a 
series  of  nodes  with  identical  excitation  conditions 
R1  = the  specified  temperature:  R1=0.0  default 
R2  « the  specified  heat  flow  rate:  R2=0.0  default 
R3  = the  specified  convective  fluid  temperature:  R3=0.0  default. 
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2)  Harmonic  Excitation/Response 


The  response  of  the  system  to  steady  harmonic  excitation  may  be  computed  by  including 
the  following  lines  of  data: 

HARMONIC 

F-R1 

11,12,13  T=R2,R3  Q=R4,R5  TEX=R6,R7 

END 

where: 

R1  = the  frequency  of  the  harmonic  excitation  [=]  cycles/time 

11,12,13  = the  first  node  number,  last  node  number,  and  node  number  increment  of  a 
series  of  nodes  with  identical  excitation  conditions 
Harmonic  Temperature  Excitation: 

R2  = Amplitude  of  Excitation:  R2  = 0.0  default 
R3  = Time  Lag:  R3  = 0.0  default 
Harmonic  Direct  (External)  Flow  Rate  Excitation: 

R4  = Amplitude  of  Excitation:  R4  = 0.0  default 
R5  = Time  Lag:  R5  = 0.0  default 

Harmonic  External  (Convective  Fluid)  Temperature  Excitation: 

R6  = Amplitude  of  Excitation:  R6  = 0.0  default 
R7  = Time  Lag:  R7  = 0.0  default 

note;  Frequency  = (circular  frequency)/27t 

Lag  is  defined  in  terms  of  the  phase  angle  or  argument  of  the  complex 
representations  of  both  excitation  components  and  response  as: 

Lag  = (-argument)^  2 n (frequency)  ) 

(See  discussion  in  section  2.6.2  of  the  REPORT,  above,  for  clarification.) 


3)  General  Excitation/Response  • Predictor-Corrector  Integration 

The  response  of  the  system,  including  transients,  to  a general  dynamic  excitation  that  is 
defined  in  terms  of  piecewise  linear  functions  may  be  computed  by  including  the  following 
lines  of  data: 

Solution  Control  Data: 


PREDICTOR-CORRECTOR 
TIME=R1,R2,R3  ALPHA=R4  NEFNsll 

12,13,14  T=I5  Q=I6  TEX=I7 
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END 


where: 


R1.R2.R3 

R4 


II 

12,13,14 

15 

16 
17 


= beginning  time,  ending  time,  time  step  increment  for  dynamic  analysis 
= integration  parameter;  0.0  < R4  < 1.0  ; R4=0.75  default;  instability  may 
result  for  R4  < 0.5  . 

= number  of  excitation  functions  that  will  be  defined. 

= the  first  node  number,  last  node  number,  and  node  number  increment  of  a 
series  of  nodes  with  identical  excitation  conditions 
= the  excitation  function  number  for  specified  temperature 
= the  excitation  function  number  for  specified  external  flow  rates 
= the  excitation  function  number  for  specified  external  fluid  temperature  for 
convective  boundary  conditions 


Initial  Conditions  Data: 


INITIAL  CONDITIONS 
11,12,13  TrRI 

END 


INITIAL  CONDITIONS 

or  RE  AD  = C1 

END 


where: 

11,12,13  = the  first  node,  last  node,  node  number  increment  of  a series  of  nodes  with 
identical  initial  temperature  conditions 
R1  = the  initial  temperature  condition:  R1=0.0  default 

Cl  = the  filename  of  file  to  read  initial  conditions  from. 


Initial  temperature  conditions  are  specified  by  as  many  lines  as  needed.  If  this  data  group  is 
omitted  initial  conditions  of  temperature  at  all  nodes  will  be  assumed  to  be  zero. 

If  the  READ=C1  option  is  used  instead,  initial  conditions  will  be  read  from  the  last  record 
of  <filename>.TMP  in  the  form  (TIME,  (T(N),  N=1,NNOD)  ) where  NNOD  is  the  number  of 
nodes  in  the  system  and  the  value  TIME  is  ignored  (see  REPORT  below).  This  option 
facilitates  restarting  a problem  to  continue  computation  or  as  a means  to  get  an  estimate  of 
initial  conditions  that  are  not  known  apriori. 

Excitation  Function  Data: 

EXCITATION 

T=R1  EFN=R2,R3,R4,  ... 


END 

where: 

R1  =time 

R1,R2,R3,  ...  = excitation  function  values  for  time=R1  for  excitation  functions  1,  2,  3, 
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...  respectively 


H.  Output  Report  Control 

The  following  lines  of  data  control  output  reporting: 

REPORT 
[WRITE]  PINTsll 
12,13,14 

END 

where: 

If  the  keyword  WRITE  is  included  the  time  and  computed  temperature  solution  for  all 
nodes  will  be  written  to  a binary  file  <filename>.TMP  as  a series  of  records  of  the  form 
(TIME,  (T(N),N=1,NNOD)),  where  NNOD  is  the  number  of  nodes  in  the  system 
II  = the  temperature  solution  "print"  output  interval  (the  solution  will  be  written  to 
a printable  ASCII  file  <filename>.OUT;  used  by  predictor-corrector  integration 
procedure  only  - otherwise  ignored 

12,13,14  = the  first  node  number,  last  node  number,  and  node  number  increment  of  a 
series  of  nodes  selected  for  printable  output. 

The  WRITE  option  is  useful,  when  used  with  the  READ  option  used  to  specify  initial 
conditions  (above)  to  restart  analysis. 
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Appendix  A - Free-Fleld  Input  Conventions  of  the  CAL-SAP  Software 
Development  System 


A "C"  in  column  1 of  any  line  will  cause  the  line  to  be  echoed  as  a comment  on  the  console. 

A backslash  V at  the  end  of  information  on  the  first  line  will  allow  the  next  line  to  be 
interpreted  as  a continuation  of  the  first  line:  therefore,  a 160  character  record  is  possible. 

A colon  indicates  the  end  of  information  on  any  line.  Information  entered  to  the  right  of 
the  colon  is  ignored  by  the  program.  Therefore,  it  can  be  used  to  annotate  an  input  file. 

Data  may  be  identified  by  an  identifier  of  the  form  "<identifier="  (e.g.  NUM=).  Data  not 
identified  by  an  identifier  must  be  placed  first  in  a line.  If  fewer  data  items  are  supplied  than 
required  the  missing  data  will  be  assumed  to  be  zero  or  blank  as  appropriate. 

Decimal  points  are  not  needed  for  real  data.  Scientific  notation  formats  of  the  form  E±  nn 
(e.g.  5.5  E-15)  may  be  used.  Simple  arithmetic  expressions  may  be  used  using  the  standard 
operators  +,  *,  *,  and  / . The  order  of  evaluation  is  sequential  from  left  to  right  - unlike 
FORTRAN. 
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Appendix  B - Example  Problems 


Bl  Response  of  a Seml-Inflnlte  Slab  to  a Step  Change  of  Surface  Heat  Flux 

To  demonstrate  the  use  and  characteristics  of  the  one-  and  two-dimensional  isoparametric 
conduction  finite  elements  the  temperature  response  of  a semi-infinite  slab,  T(x,t),  to  a step 
change  of  surface  flux,  Q(0,t)  = 10.0  Btu/sec  ; for  t > 0 , was  considered.  All  material 
properties  were  assigned  unit  values; 


Four  different  models  of  this  problem  were  considered,  each  representing  a 6 inch  deep 
solid  with  a perfectly  insulated  boundary  at  the  6 inch  depth.  These  models  are  illustrated 
below  in  Figure  Bl.  It  may  be  seen  that  these  models  are  similar  in  that  an  even  spacial 
discretization  at  0.5"  intervals  is  used  in  all  models.  Variants  of  each  of  the  three  models 
employing  one-dimensional  isoparametric  element  assemblages  based  on  different  capacity 
modeling  strategies  were  also  investigated  bringing  the  total  number  of  models  studied  to 
eight; 

Model  A1:  12  Two-node  ID  isoparametric  elements  with  "consistent"  (i.e. , exact 

evaluation  of  the  element  capacitance  matrices  using  two-point  Gauss  quadrature), 

Model  A2;  12  Two-node  ID  isoparametric  elements  with  "lumped"  element  capacitance 

matrices  (i.e.,  using  two-point  quadrature  with  quadrature  points  of  (-1.0,  +1.0)), 

Model  A3:  12  Two-node  ID  isoparametric  elements  with  "optimal"  evaluation  of  the 

element  capacitance  matrices  (i.e.,  using  two-point  quadrature  with  quadrature  points 
of  (-  V273,  + V2/3)), 

Model  Bl:  6 Three-node  ID  isoparametric  elements  with  "consistent"  (i.e.,  exact) 

1 Carslaw,  H.S.  & Jaeger,  J.C.,  Conduction  of  Heat  in  Solids.  2nd  Ed.  Oxford  University  Press, 


Conductivity  = k =1.0  Btu/sec-in-°F 
Specific  Heat  = c = 1.0  Btu-in/sec2 *-lb-°F 
Mass  Density  = p =1.0  sec2-lb/in4 


and  initial  conditions  were  assumed  to  be  zero,  T(x,0)  = 0.0  . 


The  exact  solution  to  this  problem  is  given  by1 ; 


where; 


k = k/pc 


1969. 
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evaluation  of  the  element  capacitance  matrices  (i.e.,  using  three-point  Gauss 
quadrature), 

Model  B2:  6 Three-node  ID  isoparametric  elements  with  "lumped"  element  capacitance 

matrices  (i.e.,  using  two-point  quadrature  with  quadrature  points  of  (-1.0,  0.0,  +1.0)), 

Model  Cl:  4 Four-node  ID  isoparametric  elements  with  "consistent"  (i.e.,  exact) 

evaluation  of  the  element  capacitance  matrices  (i.e.,  using  four-point  Gauss 
quadrature), 

Model  C2:  4 Four-node  ID  isoparametric  elements  with  "lumped"  element  capacitance 

matrices  (i.e.,  using  two-point  quadrature  with  quadrature  points  of  (-1.0,  -1/3,  +1/3, 

+1-0)), 

Model  D:  12  Four-node  2D  isoparametric  elements  with  "consistent"  (i.e.,  exact 

evaluation  of  the  element  capacitance  matrices  using  four-point  Gauss  quadrature). 


B - 2 


Temperature  distributions  were  computed  for  t=0.02  sec  and  t=0.10  sec.  The  results  of 
Model  A1  and  Model  D , the  consistent  ID  and  2D  isoparametric  elements  repectively,  were 
practically  identical  in  this  case  of  1 D heat  transfer.  This  is  to  be  expected  since  identical 
shape  functions  and  numerical  integration  schemes  were  used  for  these  cases.  Results  of 
analysis  for  the  seven  one-dimensional  element  assemblages  are  plotted  below. 


X (in) 


Fia.  B2  Temperature  Profiles  After  a Long  Time  Interval 

The  results  of  all  modelings  are  very  close  after  long  time  intervals,  converging  to  the  exact 
solution.  It  may  be  seen,  however,  that  the  higher  order  elements  (i.e.,  the  four-node  and 
three-node  elements)  provide  marginally  better  approximations  of  the  temperature  profiles. 
Although  it  appears  that  the  use  of  "consistent"  capacitance  modeling  for  these  higher  order 
elements  provides  slightly  better  approximations  of  the  temperature  profiles  than  the 
"lumped"  approximations  the  difference  is  slight  and  the  "lumped"  capacitance  modeling 
results  in  a diagonal  system  capacitance  matrix,  a fact  that  can  be  used  to  advantage  in 
transient  solution  algorithms. 


For  short  time  intervals  the  two-node  element  assemblages,  Models  A1-A3,  tend  to  predict 
results  that  overshoot  the  exact  solution,  as  may  be  seen  in  Fig.  B3.  Comparing  the  results 
of  the  three  capacitance  modeling  alternatives  reveals  that  the  "consistent"  capacitance 
modeling  provides  the  best  estimation  of  surface  temperature  but  overshoots  at  the  second 
node  the  most,  the  "optimal"  capacitance  modeling  provides  minimal  overall  error  in  the 
estimation  of  temperature  but  reveals  a slight  overshoot  at  the  second  node,  and  the 
"lumped"  capacitance  modeling  suffers  no  overshoot  at  the  second  node  but  proves  least 
accurate  in  modeling  the  surface  temperature.  The  lumped  capacitance  modeling  results  in  a 
model  that  is  mathematically  equivalent  to  a finite  difference  model  and  as  such  provides  a 
useful  comparison  between  finite-element  and  finite-difference  analysis  of  this  particular 
problem. 
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Fia.  B3  Temperature  Profiles  After  a Short  Time  Interval  for  Models  A1-A3 


Fta.  B4  Temperature  Profiles  After  a Short  Time  Interval  for  Model  A and  Models  B & C 

Finally,  Fig.  B4  provides  a comparison  of  the  "optimal"  two-node  element  assemblage  with 
the  three-node  assemblages,  Models  B1  & B2,  and  the  four-node  assemblages,  Models  Cl  & 
C2.  The  higher  order  elements,  again,  provide  better  estimates  of  surface  temperature  than 
the  lower  order  elements  with  the  "consistent"  higher  order  elements  out-performing  the 
"lumped"  higher  order  elements.  The  "lumped"  higher  order  elements.  Models  B2  and  C2,  on 
the  other  hand,  tend  to  exhibit  less  overshoot  and  again  lead  to  diagonal  system 
capacitance  matrices  that  can  result  in  substantial  savings  in  memory  and  solution  time. 

The  command/data  input  file  for  Model  B,  using  Gauss  quadrature  for  capacitance  modeling 
is  listed  below.  The  output  results  file  for  this  problem  is  also  listed. 
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Model  B-Gauss  Command/Data  Input  File 


DTAMl : 

Semi  Infinite  Slab:  MODEL  B-Gauss 
NODE=l 3 BAND =3 : 

ISOID : 

1 1=1, 3, 2 X=0 . 0, 1 . 0, 0 . 5 K=1 . 0 P=1 . 0 C=1 . 0 A=1.0  R=GAUSS : 

6 1=11,13,12  GEN=2  X=0 . 0, 1.0, 0.5  K=1 . 0 P=1 . 0 C=1 .0  A=1.0  R=GAUSS: 
END: 

PREDICT: 

TIME=0, 0.1, 0.005  ALPHA=0 . 75  NEFN=1 : 

1 Q=1 : 

END: 

EXCITATION: 

T=0  EFN=10 . 0 : 

T=1 . 0 EFN=10 . 0 : 

END: 

REPORT: 

PINT=4 : 
lr  7,1 
END: 

Mpdel...B-.G3uss  Output  Results  File 


| DTAMl:  Building  Energy  Simulation  for  Linear  Systems  I 

Ver-6-86-- 

Jim  Axley  - Cornell  & NBS 
MTOT : 50000 

====  PROBLEM  LABEL 
SEMI  INFINITE  SLAB 
=====  PROBLEM  CONTROL  VARIABLES 


Number  of  nodes  13 

Maximum  probable  bandwidth  ...  3 


====  EQUALITY  CONSTRAINT  CONDITIONS 

--  NOTE:  Constraint  condition  data  not  found. 

No  DOFs  will  be  constrained  to  be  equal. 

====  NODAL  COORDINATES 

--  NOTE:  Coordinate  data  not  found. 

====  element  data 

==  ISOPARAMETRIC  ID  2 TO  4-NODE  CONDUCTION  ELEMENTS 

Kx  = Conductivity  P = Density  C = Capacity 

A = Area  Q = Int.  Heat 
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Elem 

Node 

Coords 

Kx 

P 

C A 

Q' 

Quad.  Points 

1 

1 

.00 

1.0 

1.0 

1.0  1.0 

.00 

GAUSS 

3 

1.0 

2 

.50 

2 

3 

.00 

1.0 

1.0 

1.0  1.0 

.00 

GAUSS 

5 

1.0 

A 

.50 

3 

5 

.00 

1.0 

1.0 

1.0  1.0 

.00 

GAUSS 

7 

1.0 

:> 

6 

.50 

A 

7 

.00 

1.0 

1 . 0 

1.0  1.0 

.00 

GAUSS 

9 

1.0 

8 

.50 

5 

9 

.00 

1.0 

1.0 

1.0  1.0 

.00 

GAUSS 

11 

1.0 

10 

.50 

6 

11 

.00 

1.0 

1.0 

1.0  1.0 

.00 

GAUSS 

13 

1.0 

12 

.50 

— 

NOTE 

: Needed  bandwidth  

3 

Requested  bandwidth  . . 

3 

— 

NOTE 

: Time 

to  form 

equations : 

4 seconds . 

==== 

BOUNDARY . CONDITIONS 

AND  EQUATION  NUMBERS 

— 

NOTE 

: Boundary  condition  data 

not  found. 

All 

nodes  assumed  to  be 

flux-prescribed 

nodes . 

Equation  numbers  = node  numbers. 

=====  CONVECTIVE  BOUNDARY  CONDITIONS 

- — NOTE:  Convection  boundary  data  not  found. 

====  SOLUTION 

==  DYNAMIC  SOLUTION:  PREDICTOR-CORRECTOR  INTEGRATION 

--  SOLUTION  CONTROL  INFORMATION 


Start  time  .000 

End  time  .100 

Time  step  increment  0.500E-02 


Integration  parameter,  alpha  .750 

Number  of  excitation  functions  ...  1 

— EXCITATION  FUNCTION  NUMBERS 

Node  Temp  Ext. Flux  Ext. Temp. 

1 10 

— INITIAL  CONDITIONS 

• — NOTE:  Initial  condition  data  not  found. 
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Initial  temperatures  assumed  to  be  zero. 

— TIME  STEP 

--  NOTE:  Estimated  time  step  to  limit  error  to  approx.  5.00%  is: 

Specified  time  step  is:  0.500E-02 


— RESPONSE 


Time:  0.200E-01 

Node  Temp.  Node 

1 1.27  2 

6 -0.101E-03  7 

Time:  0.400E-01 

Node  Temp.  Node 

1 2.09  2 

6 0 . 142E-02  7 


Temp.  Node  Temp. 
-0 . 7 60E-01  3 .117 
0.48  6E-03 


Temp.  Node  Temp. 
0 . 1 67E-01  3 .102 

-0 . 863E-03 


Node  Temp.  Node 
4 -0.51 6E-02  5 


Node  Temp.  Node 
4 0 . 905E-02  5 


Time:  0.600E-01 

Node  Temp.  Node  Temp.  Node 

1 2.68  2 .178  3 

6 0.126E-02  7 -0.641E-03 


Temp.  Node 
0 . 678E-01  4 


Temp.  Node 
0 . 188E-01  5 


Time:  0.800E-01 

Node  Temp.  Node 
1 3.16  2 

6 0 . 4 12E-03  7 


Temp.  Node  Temp.  Node 
. 366  3 0 . 472E-01  4 

0 . 101E-03 


Temp.  Node 
0.221E-01  5 


Time:  .100 

Node  Temp.  Node 
1 3.56  2 

6 -0.21 IE-03  7 


Temp.  Node  Temp.  Node 
.564  3 0.47  6E-01  4 

0. 611E-03 


Temp.  Node 
0.216E-01  5 


--  NOTE:  Solution  time:  9 seconds. 


. 277E-02 


Temp . 

0. 932E-02 


Temp . 

-0 . 131E-03 


Temp. 

-0.44  7E-02 


Temp . 

-0. 313E-02 


Temp . 

0 . 853E-03 
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B2  Residential  Building Example 


A section  of  a hypothetical  town  house  and  a corresponding  thermal  idealization  is 
illustrated  below. 


Fig  B5  Hypothetical  Town  House  and  Thermal  Idealization 

It  may  be  seen  that  the  model  is  primarily  a R/C  network  model  with  a single  linear  one- 
dimensional two-node  isoparametric  finite  element,  LI,  used  to  model  the  dynamics  of  the 
mass  wall. 

The  command/input  file  to  compute  both  steady  and  harmonic  indoor  air  temperature 
responses  for  specified  conditions  of  outside  air  temperature  and  solar  gain  is  listed  below; 


DTAMl : 

Townhouse  Model: 

NODE=T  9 BAND=5: 

RESIST:  Consistant  Units  of  Btu,  ft,  hr 


1 

1=1,2 

A=17  R=  1/1.5 

R1 

Air  Film 

at 

Sunspace 

2 

1=2,3 

A=17 R=0 . 9 

R2 

Sunspace 

Glass 

3 

1=3,  4 

A=17  R=  1/1.5 

R3 

Air  Film 

at 

Sunspace 

4 

1=4,5 

A=8  R= 1 / 1 . 5 

R4 

Air  Film 

at 

Mass  Wall 

5 

1=6,7 

A=8  R=1 / 1 . 5 

R5 

Air  Film 

at 

Mass  Wall 

6 

1=7,  8 

A=8  R-l/1.5 

R6 

Air  Film 

at 

N.  Wall 

7 

1=8,  9 

A=8  R=ll 

R7 

N.  Wall 

8 

1=9, 10 

A=8  R-l/1.5 

R8 

Air  Film 

at 

N.  Wall 

9 

1=7,11 

A=24  R—l/1 . 1 

R9 

Air  Film 

at 

Ceiling 

10 

1=11, 12 

A=2  4 R=1 1 

RIO 

Ceiling 

11 

1=12, 13 

A=24  R—l/1 . 6 

Rll 

Air  Film 

at 

Ceiling 

12 

1=13, 14 

A=13  R=  1/1.3 

R12 

Air  Film 

at 

Roof 

13 

1=13,17 

A=13  R=l/1 . 3 

R13 

Air  Film 

at 

Roof 

14 

1=14, 15 

A=13  R=1 1 

R14 

Roof 

15 

1=17, 18 

A=13  R=1 1 

R15 

Roof 

16 

1=15, 16 

A=13  R=l/1 . 6 

R16 

Air  Film 

at 

Roof 

17 

1=18,19 

A=13  R=1 / 1 . 6 

R17 

Air  Film 

at 

Roof 

END : 

CAPACI : Consistant  Units  of  Btu,  ft,  hr 
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1 1=4  M=150*12/3  C=0.16  :C1  Sunspace  Quickmass  (4"  cone,  floor) 

2 1=7  M=150*24/3  C=0.16  :C1  Room  Quickmass  (4"  cone,  floor) 

END: 

IS01D:  Consistant  Units  of  Btu,  ft,  hr 

1 1=5,6  X=0. 0,8.0/12.0  K=1 . 0 P=150  C=0.2  A=8  :L1  8"  Concrete  Wall 

END: 


BOUNDARY : 
1 BC=T 
10  BC=T 
16  BC=T 
19  BC=T 
END: 


:Outside  air  temp  to  be  specified 
:Outside  air  temp  to  be  specified 
:Outside  air  temp  to  be  specified 
:Outside  air  temp  to  be  specified 


STEADY: 

1 T=7 . 5 

10  T=7 . 5 

16  T=7 . 5 

19  T=7 . 5 

4 Q=200*13/3.1416 

END 

HARMONIC: 


: Out side  air  temp  - Steady  Component 
: Outside  air  temp  - Steady  Component 
: Outside  air  temp  - Steady  Component 
: Outside  air  temp  - Steady  Component 
: Sunspace  Solar  Gain  - Steady  Component 


F=1 .0/24 
1 T=10 , 10 

10  T=10, 10 

16  T=10, 10 

19  T=1 0 , 1 0 

4 Q=200*13*0.5,6 

END 

REPORT : 

P INT=1 : 

If  19, 1 


: Outside  air  temp  - Harmonic  Component 
: Outside  air  temp  - Harmonic  Component 
: Outside  air  temp  - Harmonic  Component 
: Outside  air  temp  - Harmonic  Component 
: Sunspace  Solar  Gain 


END: 
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B3  Analysis  of  the  Huron  Building  Thermal  Bridge 


In  this  example,  assemblages  of  two-dimensional  isoparametric  conduction  elements  and 
simple  resistance  elements  are  used  to  model  the  details  of  heat  transfer  in  an  exterior 
masonry  wall  section,  at  the  intersection  of  a concrete  floor  slab,  Figure  B6,  of  a recently 

constructed  federal  office  building  in  Huron,  South  Dakota,  studied  earlier  by  Grot  et.  al 2 * *. 


4"  Face  Brick 
6"  Lt.  Wt.  Cone.  Block 
3"  Fiber  Glass  Insul. 
5/8"  Gypsum  Board 


Expansion  Joint 

Reinforced  Concrete 
Slab 


A A A ,A  A AAA 


■ 


Steel  Beam 


Fireproofing 
! V/7? v Air  Space 


Fia.  B6  Huron  Building  Wall  Section 


2 Grot,  R.A.,  Childs,  K.W.,  Fang,  J.B.,  & Courville,  G.E.,  "The  Measurement  and  Quantification  of 

Thermal  Bridges  in  Four  Office  Buildings",  ASHFtAE  Trans.  1985,  V.  91,  Pt.  1,  CH-85-11  No.  4 
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Field  thermographic  measurements  of  this  building  revealed  a thermal  bridge  (i.e.,  highly 
conductive  path  through  the  wall  construction)  due,  evidently,  to  the  penetration  of  the 
floor  slab  that  acted  to  significantly  increase  building  heat  loss.  Analysis  of  heat  transfer  in 
the  vicinty  of  this  thermal  bridge  is  complicated  by  the  complex  two-dimensional  nature  of 
the  geometry  and  the  air  space  between  the  steel  beam  and  the  outer  wall. 

Two  models  of  the  wall  section  were  considered.  In  the  first  model,  Model  1,  an 
assemblage  of  two-dimensional  isoparametric  elements  was  used  to  model  the  simpler 
construction  from  the  center  of  the  concrete  slab  upward,  Fig.  B8;  heat  transfer  was 
assumed  to  be  symmetric  about  the  slab  centerline.  In  the  second  model,  Model  2,  an 
assemblage  of  two-dimensional  isoparametric  elements  were  used  to  model,  in  detail,  the 
geometry  and  discontinuous  material  properties  of  the  actual  construction  and  a subassembly 
of  simple  resistance  elements  was  used  to  model  convective/radiative  heat  transfer  in  the  air 
space  adjacent  to  the  steel  beam,  Fig.  B9.  The  response  of  both  models  to  a steady  state 
temperature  drop  of  30.5°F,  from  inside-air  to  outside-air,  and  a steady  excitation  of  a 
harmonic  variation  of  outside  air  of  a 12.5°F  amplitude  and  24  hour  period  were  computed. 

The  response  of  both  models  was  similar.  Representative  results,  for  the  temperature  field 
through  the  construction,  are  shown  in  Fig.  B7.  It  is  seen  that  the  steady  state  response 
dominates  the  behavior  of  the  construction  and  the  thermal  path  provided  by  the  slab  is 
evident.  Surface  fluxes,  calculated  from  the  computed  results  closely  agreed  with  measured 
behavior  and  other  analytical  results. 


0 1 2 X 3 


Fig,  B7  Computed  Temperature.fjfild  far  Model  .1 

Model  1 employed  34  elements  connected  at  46  nodes  and  Model  2 employed  105 
elements  connected  at  139  nodes.  Computation  to  form  and  assemble  the  element 
equations  for  Model  2 took  73  seconds  and  solution  for  both  the  steady  state  and  steady 
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harmonic  excitation  took  an  additional  48  seconds  using  a Macintosh  microcomputer  with  a 
MC  68000  microprocessor  and  51 2K  of  central  memory. 
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**  Expansion  joint  assumed  to  be  equivalent  to  mineral 
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Fig,  B8  Huron  Building  Wall  .Seslign  MsdeLL 
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Appendix  C - Program  Structure 


The  program  DTAM1  is  structured  in  a manner  similar  to  a command  processor. 
Structurally  it  was  developed  as  a collection  of  trees  of  subroutines  linked  to  the  main 
program  through  roots . 


Fia.  Cl  Generic  Trees  Structure 


Each  tree  was  developed  to  be  largely  independent,  dynamically  defining  arrays  and 
seeking  only  the  input  data  that  it  alone  requires. 


Root  of  Tree:  - dynamically  DEFINE's  arrays  used  by  Tree 


Tree:  - gets  only  input  data  needed  by  tree 
- executes  command  operation 

Locally  shared  labeled  COMMON  & SUBROUTINES 


Fig.  C2  Generic  Tree  of  Subroutines 

This  independence  allows  incremental  processing  of  input  data  and  reporting  of  results  so 
that  a single  error  in  the  command/data  input  data  file  will  not  result  in  program  termination. 
It  is  also  hoped  that  the  hierarchical  structure  will  help  to  facilitate  future  program 
development  efforts. 

The  CAL-SAP  development  software  subroutines  [2]  are  used  by  all  trees  for  dynamic  array 
management  and  input  data  interpretation.  The  CAL-SAP  dynamic  array  management  is 
accomplished  by  storing  array  values  and  a directory  to  those  values  in  a single  blank 
common  vector  IA(MTOT)  as  indicated  below  in  Fig.  C3.  One  may  change  the  array  capacity 
of  the  program  by  simply  resetting  the  variable  MTOT  and  redimensioning  IA(MTOT)  in  the 
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main  program  declaration. 


Arrays  stored  low 
inblark  common. 


Directory  of  arrays  with  pointers 
to  arrays  stored  high  in  blank  common. 


t* 

IA(1) 


[A]  {B} 


1C] 

(D) 

DIR  D 

DIR  C 

DIR  B 

DIR  A 


■MTOT- 


4 


IA(MTOT) 


Fia.  C3  CAL-SAP  Dynamic  Array  Management 


Specifically,  the  main  program  DTAM1  calls  six  trees  of  subroutines; 


Fia.  C4  Structure  of  DTAM1 

RID  is  the  root  to  the  ID  Tree  that  establishes  node-number  to  equation-number 
correspondance,  accounting  for  any  constraint  input;  it  creates  an  ID  array. 

RCOORD  is  the  root  to  the  Coordinate  Tree  that  inteprets  coordinate  input  data  and 
generates  nodal  coordinates;  it  forms  the  XYZ  array. 

RFORM  is  the  root  to  the  FORM  Tree  that  forms  element  arrays  and  assembles  them  to 
form  system  arrays;  it  forms  the  system  conductance  matrix,  K(NEQN.MBAND), 
the  system  capacitance  matrix,  C(NEQN,MBAND),  and  the  internal  generation 
heat  flow  rate  vector,  EO(NEQN)  - where  NEQN=the  number  of  equations  and 
MBAND=the  system  (half)  bandwidth. 

RBOUND  is  the  root  to  the  BOUND  Tree  that  establishes  boundary  conditions  by 
modifying  the  ID  array. 

RCONV  is  the  root  to  the  CONV  Tree  that  establishes  convective  boundary  conditions  by 
modifying  the  system  conductance  matrix  and  stores  convective  coefficients  in 
array  CH(NNOD),  where  NNOD  = the  number  of  nodes  in  the  system. 

RLSOLV  is  the  root  to  the  LSOLV  Tree  that  interprets  solution  specification  and 
excitation  data  and  implements  solution  by  calling  STEADY,  HARMON,  and/or 
PREDIC  to  affect  steady  state,  steady  harmonic,  or  predictor-corrector  solution 
options. 
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APPENDIX  D 

FORTRAN  77  SOURCE  CODE 


C 


The  FORTRAN77  source  code  for  DTAM1  is  listed 
below. 


PROGRAM  DTAM1  C LINEAR:  INTERNAL  VARIABLES 

ERR*. FALSE. 

C — PRO: DTAM1  - A BUILDING  ENERGY  SIMULATION  PROGRAM  BASED  ON  A C 

C DISCRETE  THERMAL  ANALYSIS  METHOD  FOR  LINEAR  IDEALIZATIONS  C— 2.0  WRITE  BANNER  & OPEN  INPUT  AND  OUTPUT  DATA  FILES 

C C 

c DEVELOPED  BY  JAMES  AXLEY . CORNELL  UNIVERSITY  WRITE (NTM,  2200)  MTOT 

C NBS  BUILDING  PHYSICS  DIV.  1985-86  CALL  IFILE 

C USING;  CALL  POPENCOUT  ’) 

C A)  CAL-SAP  LIBRARY  OF  SUBROUTINES  DEVELOPED  BY  ED  WILSON.  WRITE (NOT. 2200)  MTOT 

C U.C.  BERKELEY 

C B)  MacFort ran  V2 . 0 4 Microsoft  Fortran  V2.1  COMPILERS  2200  FORMAT (// 

C ANSI  FORTRAN  77  STANDARD  .*  

C • INCLUDING;  • ' / 

C 1.  USE  OF  SAVE  /contnonname/  TO  RETAIN  DEFINITION  OF  .*1  DTAM1 : Building  Energy  Simulation  for  Linear  Systems 

C LOCAL,  SUBROUTINE  VARIABLES  AND  SOME  NAMELIST  VARIABLES  . I'/ 

C ‘EXCEPT;  . * 

C 1.  TYPE  MacFort ran  EXTENSION  "TYPE"  USED  FOR  . Ver-6-86— •/ 

C PROMPTS  TO  ALLOW  IN-LINE  RESPONSES  .55X,’Jlm  Ax ley  - Cornell  4 NBS*/ 

C*« •••*••*•••••••••••*•••  . 65X, 'MTOT: ’110) 

c c 

IMPLICIT  REAL* 8 (A-H.O-Z)  C--3.0  WHILE  ERR*. FALSE. 

C C 

C****  CAL-SAP:  DATA  4 COMMON  STORAGE  **•*•*•*•*••*•**••*•****•••••••**  C 3.1  GET  PROB  LABEL  4 CONTROL  INFORMATION 

C C 


CAL-SAP:  LOGICAL  UNIT  NUMBERS 

NTM  - 9 
NTR  - 9 
NIN  « 10 
NOT  = 11 
ND1  « 12 
ND2  - 13 
ND3  - 14 
ND4  * 15 


CHARACTER  FIN*1.EXT*1 
COMMON  MTOT. NP. IA (50000) 

COMMON  /DBSYS/  NUMA . NEXT. IDIR, IP (3 ) 

COMMON  /IOLIST / NTM, NTR, NIN, NOT, ND1 , ND2 , ND3 . ND4 
COMMON  /PARC/  FIN  (12) . EXT (8  0) 

C 

C DICTIONARY  OF  VARIABLES 

C 

C VARIABLE 

C MTOT 

C NP 

C IA (MTOT) 

C NUMA 

C NEXT 

C IDIR 

C IP  (3) 

C NTM 

C NTR 

C NIN 

C NOT 

C ND1  thru  ND4 

C FIN  (12) 

C EXT (80) 

c**..»**»*t..t*t*«i 

C 

C****  LINEAR:  DATA  4 COMMON  STORAGE 
C 


DESCRIPTION 

SIZE  OF  BLANK  COMMON  VECTOR  IA 

CURRENT  DATA  TYPE:  1 "INTEGER;  2= REAL ; 3*LOGICAL 
BLANK  COMMON  VECTOR 

NUMBER  OF  ARRAYS  IN  BLANK  COMMON  DATA  BASE 
NEXT  AVAILABLE  STORAGE  LOCATION  IN  BLANK  COMMON 
START  OF  DIRECTORY  IN  BLANK  COf**ON 
NUMBER  OF  LOGICALS  IN  EACH  DATA  TYPE 
LOGICAL  UNIT  NUMBER  FOR  TERMINAL /SCREEN 
LOGICAL  UNIT  NUMBER  FOR  TERM INAL /KEY BOARD 
LOGICAL  UNIT  NUMBER  FOR  INPUT  DATA  FILE 
LOGICAL  UNIT  NUMBER  FOR  OUTPUT  DATA  FILE 
LOGICAL  UNIT  NUMBERS  FOR  GENERAL  USE 
INPUT  DATA  FILE  NAME 

FILE  EXTENSION  NAME  FOR  GENERAL  USE 


LOGICAL*l 

CHARACTER 


ERR 

YESNO*l . SAVE  *4 


C DICTIONARY  OF  VARIA 


L E S 


VARIABLE  DESCRIPTION 

ERR  DO-WHILE  TERMINATOR  FLAG 

YESNO  CONTINUATION  FLAG 

NNOD  NUMBER  OF  NODES  IN  SYSTEM 

NEON  NUMBER  OF  (SYSTEM)  EQUATIONS 

MBAN  (HALF)  BANDWIDTH  OF  SYSTEM  EQUATIONS 

POINTERS  TO  BLANK  COMMON  LOCATIONS 


C MPXYZ 

C MPC 

C MPK 

C MPEO 

VECTOR 
C MPT 

C MPID 

ARRAY 
C 
C 

NODE 

C 

C MPCH 

FLUX 

C MPNPR 

C 


XYZ (NNOD , 3)  : COORDINATE  ARRAY 

C(NEQN.MBAN) : CAPACITY  MATRIX 

K (NEON. MBAN)  : CONDUCTANCE  MATRIX 

EO  (NNOD)  : INTERNALLY  GENERATED  EXCITATION 


T (NEQN) 
ID (NNOD) 


CH  (NNOD) 
NPR  (NNOD) 


NODAL  TEMPERATURE  VECTOR 

NODAL  EQUATION  NUMBERS  B.  C.  FLAG 


EQTN  NO.  « ABS (ID) 

WEG.  ID  - TEMPERATURE  PRESCRIBED 


POS.  ID  - FLUX  PRESCRIBED  NODE 
COEFS  TO  COMPUTE  EFFECT.  CONVECT. 


OUTPUT  PRINT  CONTROL  ARRAY 


C — 1.0  INITIALIZE  INTERNAL  VARIABLES 


C CAL-SAP:  DATABASE 

10  MTOT  * 50000 
NUMA  - 0 
NEXT  * 1 
IDIR  = MTOT 
IP(1)  ■ 4 
IP  (2)  * 8 
IPO)  « 1 


CALL  FINDN ( ' DTAM1 ’ , 5 . KEY ) 

IF  (KEY  . EQ.  1 ) THEN 
WRITE (NTM. 2310) 

WRITE (NOT, 2310) 

GO  TO  999 

ENDIF 

2310  FORMAT ( ' *•*•  ERROR:  Separator  "DTAMl * not.  found.’) 

CALL  FREETY 
WRITE (NTM. 2312) 

WRITE (NOT. 2312) 

2312  FORMAT  (/’  **=  PROBLEM  LABEL’/) 

CALL  FREE 
CALL  FREEPT 

WRITE (NTM. 2316) 

WRITE (NOT,  2316) 

2316  FORMAT  (/’  «“=  PROBLEM  CONTROL  VARIABLES’) 

NNOD=0 
MBAN=0 
CALL  FREE 
CALL  FREETY 

CALL  FREEI ('E’ .NNOD. 1) 

IF  (NNOD. LE. 0)  THEN 

WRITE (NTM, 2314) 

WRITE (NOT. 2314) 


ERR= . TRUE . 

ENDIF 

2314  FORMAT ( ’ •*•*  ERROR:  Number  of  nodes  must  be  greater  than  0.’) 
CALL  FREEI (’D’ . MBAN. 1) 

IF (MBAN .LE . 0)  MBAN  - 5 + INT (SQRT (REAL (NNDO) ) ) 


WRITE (NOT.  2318)  NNOD, MBAN 
2318  FORMAT (/ 

.*  Nurcber  of  nodes  *,I5./ 

.'  Maximum  probable  bandwidth  ...*,15) 


CALL  FREE 

CALL  FREEH ('  ’.SAVE, 4.1) 

IF (ERR)  OO  TO  999 
C 

C 3.2  ESTABLISH  EQUATION  NUMBERS 

C 

CALL  RID  (MPID, NNGD, NEQN, ERR) 

IF (ERR)  OO  TO  999 
C 

C 3.3  GET  4 GENERATE  NOOAL  COORDINATES 

C 

CALL  RCOORD  (MPXYZ.  NNOD,  ERR) 

IF (ERR)  OO  TO  999 
C 

C 3.4  FORM  ELEMENT  ARRAYS  AND  ADD  TO  SYSTEM  ARRAYS 

C 

CALL  TIME  (ITIME1 ) 

CALL  RFORM (MPID, MPXYZ . MPC , MPK. MPEO, NNOD. NEQN . MBAN . ERR) 

CALL  TIME (ITIME2) 

WRITE  (NTM. 2340)  (ITIME2-ITIME1 ) 

WRITE (NOT. 2340)  (ITIME2-ITIME1 ) 

2340  FORMAT(/'  — NOTE:  Time  to  form  equations:  *,I5,'  seconds.*) 


D -1 


IF (ERR)  GO  TO  999 
C 

C 3.5  GET  BOUNDARY  CONDITION  FLAGS 

C 

CALL  RBCOND (MPID. NNOD. ERR) 

IF (ERR)  GO  TO  999 
C 

C 3.6  MODIFY  K (NEQN.  MBAN)  FOR  CONVECTIVE  TRANSFER 

C 

CALL  RCONV (MPID. MPXYZ.  HPK, MPCH . NNOD, NEQN , MBAN . ERR) 

IF  (ERR)  GO  TO  999 
C 

C 3.7  SAVE  SYSTEM  ARRAY  DATABASE 

C 

IF (SAVE. EQ. ’SAVE’ ) CALL  SAVE 
C 

C 3.8  SOLVE  SYSTEM  OF  EQUATIONS 

C 

CALL  RLSOLV (MPID, MPK. MPC, MPEO. MPCH, NNOD. NEQN, MBAN, ERR) 

C 

C— 4.0  CLOSE  INPUT  AND  OUTPUT  DATA  FILES 
C 

999  CALL  FCLOSE  (NOT) 

CALL  FCLOSE  (NIN) 

C 

C — 5.0  SOLICIT  ANOTHER  PROBLEM 
C 

TYPE  2500 

READ (NTM, 2510)  YESNO 

IF ( (YESNO.EQ. *Y') .OR. (YESNO.EQ. *y’ ) ) GO  TO  10 
2500  FORMAT (/'  ••  Do  you  want  to  consider  another  problem?  (Y/N)  ’) 

2510  FORMAT (LAI) 

STOP 

END 

C*M*t****t**..M«t****.t«t****t*tt*t*«**««*l.**M*t*«t*M**Mtt*.**.* 

C TREE:  ID  - TREE  OF  SUBROUTINES  TO  ESTABLISH  EQUATION  NUMBERS 


C RID 

SUBROUTINE  RID (MPID. NNOD, NEQN. ERR) 

C— SUB: RID  - ROOT  TO  ESTABLISH  THE  NODE  NUMBER-EQUATION  NUMBER 
C CORRESPOND AN CE 

C PRESENTLY  IMPOSES  EQUALITY  CONSTRAINTS  AND  NUMBERS  EQUATIONS 

C *••  IN  NODE  NUMBER  ORDER 
C ***  TO  BE  USED  SUBSEQUENTLY  TO; 

C •**  1.  OPTIMIZE  EQUATION  NUMBERING  TO  MINIMIZE  BANDWIDTH 

C 2.  ALLOW  USER  REDEFINITION  OF  NUMBERING 

C 

C CAL-SAP:  DATA  6 COMMON  STORAGE 

C 

COMMON  MTOT.NP, IA (1) 

COMMON  /DBSYS / NUMA . NEXT . IDIR, IP  (3 ) 


COMMON  /IOLIST/  NTM, NTR,  NIN , NOT. ND1 , ND2 , ND3 . ND4 
C 

C COORD:  DATA  t COMMON  STORAGE 

C 

LOGICAL* 1 ERR 

SAVE  /DBSYS/. /IOLIST/ 

C 

C — 0.0  WRITE  HEADER 
C 

WRITE (NTM, 2000) 

WRITE (NOT. 2000) 

2000  FORMAT (/'  = EQUALITY  CONSTRAINT  CONDITIONS'/) 

C 

C— 1.0  DEFINE  ID  ARRAY  & INITIALIZE 
C 

CALL  DELETE (' ID  ') 

CALL  DEFINI ('ID  MPID . NNOD. 1 ) 

CALL  ZEROI (IA  (MPID) , NNOD. 1) 

C 

C 2.0  FIND  SEPARATOR  <CGNSTR>AIN;  NUMBER  IN  NODE  ORDER  IF  NOT  FOUND 

C 

CALL  FINDN ( 'CONSTR* , 6,  KEY) 

IF (KEY . EQ. 1)  THEN 
WRITE  (NTM,  21 00) 

WRITE  (NOT.  2100) 


C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL-SAP:  DATA  6 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR. NIN , NOT . ND1 . ND2 , ND3 . ND4 
C 

C ID:  DATA  i COMMON  STORAGE 

C 

INTEGER  MSTR,  MNODE.  MDOF 
CHARACTER  ENDFLAG*  3 , BC*1 
DIMENSION  ID  (NNOD),  IJK(3) 

LOGICAL* 1 ERR 

SAVE  /IOLIST/ 

C 

C — 1.0  READ  CONSTRAINT  DATA  AND  GENERATE  CONSTRAINT  FLAGS 
C 

10  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH  ( ' ' .ENDFLAG. 3.1) 

IF (ENDFLAG. EQ. 'END' ) GO  TO  20 

C GET  MASTER  NODE 

CALL  FREEI ('M' .MSTR, 1) 

IF (MSTR. LE.O. OR. MSTR. GT. NNOD)  THEN 
WRITE (NTM. 2100)  MSTR 

WRITE (NOT. 2100)  MSTR 

ERR  - .TRUE. 

ENDIF 

2100  FORMAT (’  ••*•  ERROR:  Master  node* ,15.*  Is  out  of  range.’) 

C GET  SLAVE  NODE  GENERATION  DATA 

CALL  FREEI ('  ’.IJK(1).3) 

IF (I JK (2) .EQ. 0)  IJK(2)-IJK(1) 

IF  (I  JK  (3)  . EQ.  0)  I JK  (3)  =1 
WRITE(NOT. 2110)  MSTR.IJK 

2110  FORMAT (’  Master  Node  - ',13,  • Slave  Nodes  = ' , 13. ' to  ',13 

♦ * Step  ’,13) 

DO  12  N-IJK(l) ,IJK(2).IJK(3) 

IF  (N . LE . 0 . OR. N . GT . NNOD)  THEN 
WRITE (NTM. 2120)  N 
WRITE (NOT, 2120)  N 
ERR  «=.  TRUE. 

GO  TO  10 
ENDIF 

12  ID  (N ) = -MSTR 
GO  TO  10 

2120  FORMAT ('  *•*•  ERROR:  (Generated)  Node '.15.'  is  out  of  range.’) 

C 

C--2.0  NUMBER  UNCONSTRAINED  DOF  (ID=0) 

C 

20  IF (ERR)  RETURN 
NEQN  *=  0 
DO  22  N*=l,NNOD 
IF (ID (N) . EQ. 0)  THEN 
NEQN  = NEQN  ♦ 1 

ID  (N)  «=  NEQN 
ENDIF 

22  CONTINUE 
C 

C--3.0  NUMBER  CONSTRAINED  DOF  (ID=NEG)  TO  MASTER  NODE  (WODE)  DOF 
(MDOF) 

C 

DO  30  N=1 , NNOD 
IF (ID (N) .LT. 0)  THEN 
MNODEeABS (ID (N) ) 

MDOF  - ID(WODE) 

IF (MDOF.LT. O)  THEN 

WRITE (NTM. 2300)  N. MNODE 

WRITE (NOT.  2300)  N,  MNODE 

ERR  « .TRUE. 


NEQN  - NNOD 

DO  20  N-MPID. MPID^NNOO 
20  IA  (N)  - M-HPID+1 
RETURN 
ENDIF 

2100  FORMAT ( * — NOTE:  Constraint  condition  data  not  found.*/ 

.*  No  DOFs  will  be  constrained  to  be  etjual.') 

CALL  FREETY 

C 

C— 3.0  CALL  ID  TO  DO  THE  WORK 
C 

CALL  ID  (IA  (MPID)  .NNOD.  NEQN.  ERR) 


ELSE 

ID  (N)  - »OF 
ENDIF 

ENDIF 

30  CONTINUE 

2300  FORMAT ( ' **•*  ERROR:  Slave  node  ',15.'  Is  constrained  to 

♦ slave  node  ',15,/ 

♦ ' Multilevel  constraints  are  not  permitted.') 

RETURN 

END 


RETURN 

END 


C ID 

SUBROUTINE  ID (ID. NNOO, NEQN . ERR) 

C-SUB : ID  - READS  AND  GENERATES  EQUALITY  CONSTRAINT  DATA 
C AND  NUMBERS  DOFS  VIA  NODE  ID  ARRAY 


C********************************************************************* 

C TREE:  COORD  - TREE  OF  SUBROUTINES  TO  FORM  XYZ  COORDINATE  ARRAY 

C****ttt«t***M*.****tM*t«*MMM»*ttt**t*M*t*tMt***t*Mt*.tt.ttt*. 

C R COORD 

SUBROUTINE  RCOORD (MPXYZ  .NNOD. ERR) 

C— SUB : RCOORD  - ROOT  SUBROUTINE  TO  COORD 

C COORD  - READS  COORDINATE  INFORMATION  AND  FORMS  X C Y ARRAYS 

C 


D -2 


c dictionary  of  variables 


C VARIABLE  DESCRIPTION 

C VARIABLES  PASSED  TO  SUBROUTINE 
C ERR  DO-WHILE  TERMINATOR  FLAG 

C NNOD  NUMBER  OF  NODES  IN  SYSTEM 

C MPXYZ  POINTER  TO  XYZ  (NNOD)  IN  BLANK  COMMON 

C XYZ  (NNOD)  COORDINATE  ARRAY  (ORDERED  BY  NODE  NUMBER) 


C- 


C- CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  MTOT.NP. IA (1) 

COMMON  /DBSYS / NUMA. NEXT. IDIR. IP (3 ) 

COMMON  /IOLIST/  NTM, NTR, NIN , NOT, ND1 . ND2 , ND3 , ND4 
C 

C COORD:  DATA  4 COMMON  STORAGE 

C 

LOGICAL*l  ERR 

SAVE  /DBSYS/, /IOLIST/ 

C 

C--0.0  WRITE  HEADER 
C 

WRITE (NTM, 2000) 

WRITE (NOT, 2000) 

2000  FORMAT (/'  «=  NODAL  COORDINATES ' ) 

C 

C— 1.0  FIND  SEPARATOR  <COOR>DINATES 
C 

CALL  FINDN  ( 'COORD  ' , 5,  KEY) 

IF  (KEY . EQ. 1)  THEN 
WRITE  (NTM.  2100) 

WRITE  (NOT.  2100) 

RETURN 


N2  - NN  (2) 

N3  - NN  (3) 

IF  (N1 . EQ . N2)  GO  TO  300 
IF  (N2  . EQ . 0 ) GO  TO  300 
IF (N2 . EQ.N1)  GO  TO  300 
IF  (N3  . EQ.  0)  N3  = 1 
IF  (N1 .GT.NNOD)  N1  = NNOD 
IF  (N2.GT.NNOD)  N2  » NNOD 
NDIF  « (N2-N1) /N3 
DIF  « DBLE(NDIF) 

XDIF  - ( XYZ (N2 , 1 ) - XYZ(Nl.l)  ) /DIF 

YDIF  « ( XYZ (N2, 2)  - XYZ(N1,2)  ) /DIF 

ZDIF  = ( XYZ (N2 , 3)  - XYZ(N1,3)  ) /DIF 

C GENERATE  ADDITIONAL  COORDINATES  

XX  « XYZ  (Nl.  1) 

YY  - XYZ  (Nl , 2) 

ZZ  ■ XYZ  (Nl , 3) 

Nl  - Nl  ♦ N3 
N2  - N2  - N3 

DO  200  J-N1.N2.N3 
XX  - XX  ♦ XDIF 
XYZ ( J, 1)  - XX 
YY  « YY  ♦ YDIF 
XYZ ( J, 2)  « YY 
ZZ  - ZZ  ♦ ZDIF 
XYZ ( J, 3)  - ZZ 
200  CONTINUE 

300  CONTINUE 

C PRINT  FINAL  COORDINATES  

400  DO  500  I-l.NNOD 
X - XYZ (1,1) 

Y - XYZ  (I,  2) 

Z « XYZ (1.3) 

WRITE (NOT. 2000)  I.X.Y.Z 
500  CONTINUE 
C 

900  RETURN 


ENDIF 

CALL  FREETY 
C 

C — 2 . 0 DEFINE  ARRAY 
C 

CALL  DELETE (’XYZ  ’) 

CALL  DEFINE  (’XYZ  MPXYZ, NNOD. 3) 

C 

C — 3 . 0 INITIALIZE  ARRAYS 
C 

CALL  ZEROR (IA (MPXYZ) .NNOO, 3) 

C 

C—  4.0  CALL  COORD  TO  DO  THE  WORK 
C 

CALL  COORD (IA  (MPXYZ) .NNOD, ERR) 

RETURN 

2100  FORMAT (/*  — NOTE:  Coordinate  data  not  found.’) 

END 

C 

COORD 

SUBROUTINE  COORD (XYZ. NNOD. ERR) 

C 

C SUBROUTINE  TO  READ  AND  GENERATE  COORDINATES 
C DEVELOPED  BY  MARC  HOIT  - U.C.  BERKELEY 

C 

C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

CHARACTER  END FLAG *3 
DIMENSION  XYZ (NNOD. 3) .NN (3) 

COMMON  /IOLIST/  NTM, NTR , NIN . NOT, ND1 , ND2 . ND3 . ND4 
SAVE  /IOLIST/ 

C 

WRITE (NOT, 1000) 

C INITIALIZATION 

X - 0.0 
Y - 0.0 
Z = 0.0 
NODE  - 0 

C READ  LINE  OF  COORDINATE  INFORMATION 

DO  300  I-l.NNOD 
CALL  FREE 

CALL  FREETY 

CALL  FREEH (’  ’ , ENDFLAG, 3 , 1) 


c FORMAT  SPECIFICATIONS  

C 

1000  FORMAT (/ 

1 SX, ’ NODE  X -COORD.  Y -COORD . Z -COORD. ’ ) 

C 


2000  FORMAT (7X, 13. 9X.G10.3. 4X, G1 0 . 3. 4X. G10 . 3) 
C 

END 


TREE:  FORM  - TREE  OF  SUBROUTINES  TO  FORM  AND  ASSEMBLE  ELEMENT  EQTNS 

RFORM 

SUBROUTINE  RFORM (MP ID. MPXYZ . MPC, MPK, MPEO, NNOD , NEQN. MBAN. ERR) 
--SUB:  RFORM  - ROOT  SUBROUTINE  TO  ALL  ELEMENT  SUBOUTINES  THAT 
FORM  4 ASSEMBLES  ELEMENT  ARRAYS 

DICTIONARY  OF  VARIABLES 


VARIABLE 
VARIABLES  PASSED 
ERR 
NNOD 
NEON 
MBAN 
MPEO 
MPK 
MPC 
MPXYZ 
MPID 
EO  (NNOD) 

C (NEQN. MBAN) 
K(NEQN. MBAN) 
ID (NNOD) 


XYZ (NNOD. 3) 


DESCRIPTION 

TO  SUBROUTINE 
DO-WHILE  TERMINATOR  FLAG 
NUMBER  OF  NOOES  IN  SYSTEM 
NUMBER  OF  (SYSTEM)  EQUATIONS 
(HALF)  BANDWIDTH  OF  SYSTEM  EQUATIONS 
POINTER  TO  EO (NNOD)  IN  BLANK  COMMON 
POINTER  TO  K (NEQN. MBAN)  IN  BLANK  COMMON 
POINTER  TO  C (NEQN. MBAN)  IN  BLANK  COMMON 
POINTER  TO  XYZ (NNOD. 3)  IN  BLANK  COMMON 
POINTER  TO  ID (NNOD)  IN  BLANK  COMMON 
INITIAL  EXCITATION  VECTOR 
SYSTEM  CAPACITY  MATRIX  (COMPACT  FORM) 

SYSTEM  CONDUCTANCE  TRANSFER  MATRIX  (COMPACT  FORM) 
EQUATION  NUMBER/B.C.  FLAG  ARRAY; 

ABS  (ID  (N) ) = EQUATION  NUMBER  OF  NODE  N 
SIGN (ID  (N) ) - -1  - TEMPERATURE  PRESCRIBED  NODE 
SIGN (ID (N))  - +1  - FLUX  PRESCRIBED  NODE 
COORDINATE  ARRAY  (ORDERED  BY  NODE  NUMBER) 


CAL-SAP:  DATA  4 COMMON  STORAGE 


COMMON  MTOT.NP, IA(1) 

COMMON  /DBSYS/  NUMA. NEXT. IDIR. IP (3) 

COMMON  /IOLIST/  NTM. NTR. NIN . NOT. ND1 . ND2 . ND3 . ND4 
C 

C RFORM:  DATA  4 COf*)ON  STORAGE 

C 

LOGICAL*!  ERR.  NOELEM 


IF (ENDFLAG. EQ. ’END’)  OO  TO  400 
CALL  FREEK’  ’.NODE.l) 

IF (NODE.EQ. 0)  GO  TO  400 
CALL  FREER  ( »X\ X,  1) 

CALL  FREER (’Y’.Y.l) 

CALL  FREER ( ' Z * , Z. 1) 

C STORE  COORDINATES 

XYZ  (NODE, 1)  - X 
XYZ  (NODE. 2)  « Y 
XYZ  (NODE.  3)  «=  Z 

C READ  AND  SET  GENERATION  PARAMETERS 

NN  (1)  - 0 
NN  (3)  -0 

CALL  FREE I ( 'N ’ , NN , 3) 

Nl  - NN  (1 ) 

IF(Nl.EQ.O)  GO  TO  300 


SAVE  /DBSYS/, /IOLIST/ 

C 

C— 0.0  WRITE  HEADER 
C 

WRITE  (NTM.  2000) 

WRITE (NOT. 2000) 

2000  FORMAT  (/’  -=  ELEMENT  DATA*) 

C 

C— 1.0  DEFINE  ARRAYS 
C 

CALL  DELETE  (’C  *) 

CALL  DELETE (’K  ') 

CALL  DELETE  CEO  ') 

CALL  DEFINE ( ’ EO  ’. MPEO, NNOD. 1) 
CALL  DEFINECK  ’.  MPK.  NEQN . MBAN ) 

CALL  DEFINE (’C  ’, MPC, NEQN . MBAN ) 


D - 


3 


c 

C— 2.0  INITIALIZE  ARRAYS 
C 

CALL  ZEROR  (IA  (KPEO)  .NNOO.  1) 

CALL  ZEROR  (IA  (MPK)  , NEON , MEAN) 

CALL  ZEROR  (IA  (KPC)  , NEQN  , MBAN ) 

C 

C— 3.0  FIND  ELEMENTS 
C 

NOE LEM  -.TRUE. 

MAXBAN  - 0 

C 

C— 3.1  FIND  SEPARATOR  "RESIST" 

C 

CALL  FINDN ( 'RESIST' , 6, KEY) 

IF (KEY . EQ. 0)  THEN 
WQRLEM- . FALSE . 

CALL  FREETY 

CALL  RES  I ST  (IA  (HP  ID)  . IA  (MPK)  . NNOD, NEQN . MBAN,  MAXBAN.  ’RES'  , ERR) 
END  IF 
C 

C — 3 . 2 FIND  SEPARATOR  "rLOWLO" 

C 

CALL  FINDN ( 'FLOWLO* , 6,  KEY) 

IF (KEY . EQ. 0)  THEN 
JKDELEN= . FALSE . 

CALL  FREETY 

CALL  RESIST  (IA(MPID)  . IA  (MPK)  , NNOD.  NEQN  . MBAN . MAXBAN . ’FLO’  .ERR) 
END  IF 

C 

C— 3.3  FIND  SEPARATOR  - CAP  AC  I " 

C 

CALL  FINDN  (’CAP  ACI  ’ .6.  KEY) 

IF(KET.EQ.O)  THEN 
■GELEM= . FALSE . 

CALL  FREETY 

CALL  CAPACI  (IA(MPID)  .IA(MPC)  . NNOD.  NEQN  . MBAN.  ERR) 

END  IF 
C 

C— 3.4  FIND  SEPARATOR  "ISOID" 

C 

CALL  FINDN  (*  ISOID’ .5,  KEY) 

IF  (KEY . EQ. 0)  THEN 
N€ELEM= . FALSE . 

CALL  FREETY 

CALL  ISOLD  (IA(MP ID)  .IA(MPXYZ)  .IA(MPC)  , IA  (MPK)  , IA (MPEO)  , 

. NNOD.  NEQN.  MBAN,  MAXBAN,  ERR) 

ENDIF 

C 

C— 3.5  FIND  SEPARATOR  "IS02D4" 

C 

CALL  FINDN  ( ’IS02D4  ' , 6.  KEY) 

IF  (KEY . EQ. 0)  THEN 
SCELEM* . FALSE . 

CALL  FREETY 

CALL  IS02D4  (IA  (MPID)  . IA  (MPXYZ)  . IA  (MPC)  , IA  (MPK)  , IA  (MPEO)  . 

. NNOO.  NEQN. MBAN.  MAXBAN.  ERR) 

ENDIF 

C 

C— 3.6  FIND  SEPARATOR  "VNMRT" 

C 

CALL  FINDN  ('VNMRT'  .5,  KEY) 

IF (KEY .EQ. 0)  THEN 
ROELEM=.  FALSE. 

CALL  FREETY 

CALL  VNMRT  (IA  (MPID)  , IA(MPK)  . NNOD.  NEQN , MBAN . MAXBAN . ERR; 

ENDIF 

C 

C — 4 . 0 IF  BO  ELEMENTS  ... 

C 

XFCNOCLEM)  THEN 
WRITE  (NTH.  2400) 

WRITE  (BOT, 2400) 

ERR-.TRDE. 

ENDIF 

C 

C — 5 . 0 REPORT  BANDWIDTH  NEEDED 
C 

WRITE (NTM. 2500)  MAXBAN . MBAN 
WRITE (NOT, 2500)  MAXBAN. MBAN 

RETURN 

2400  FORMAT ( ' ERROR:  No  element  data  found.) 

2500  FORMAT </’  — NOTE:  Needed  bandwidth ',15./ 

Requested  bandwidth  ..’,15) 

END 


RESIST 

SUBROUTINE  RESIST (ID, K. NNOD. NEQN . MBAN . MAXBAN . TYPE. ERR) 

C-SUB : RESIST  - FORMS  4 ASSEMBLES  RESISTANCE  OR  FLOWLOOP  ELEMENTS 
C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR, NIN . NOT . ND1 . ND2 , ND3 . ND4 
C 

C RESIST:  DATA  4 COW  ON  STORAGE 

C 

REAL *8  K (NEQN, MBAN) 

INTEGER  ID (NNOD) , LMNEN  (2) , LMOLD (2) 

LOGICAL* 1 ERR 

CHARACTER  ENDFLAG* 3 . TYPE*3 

SAVE  KNEW, NOLD, LMNEM, LMOLD. R. A, /IOLIST/ 

NDOF  - 2 
C 

C — 1.0  WRITE  ELEMENT  HEADER 
C 

IF (TYPE. EQ. 'RES’ ) THEN 
WRITE (NOT, 2100) 

ELSEIF (TYPE . EQ.  ' FLO' ) THEN 
WRITE (NOT. 2110) 

ENDIF 

2100  FORMAT ( 

. /’  — RESISTANCE  ELEMENTS’// 

.'  Elern  I-Node  J-Node  Resist.  Area') 

2110  FORMAT ( 

./’  -=  FLOW-LOOP  ELEMENTS’// 

.*  Elern  I-Node  J-Node  Flow  Rate  Capac.') 

C 

C — 2.0  GET  FIRST  ELEMENT,  FORM  4 ADD  ELEMENT  TO  SYS 
C 

INCR  « 0 

R - 0.0 

A - 0.0 

CALL  FREE 

CALL  FREETY 

CALL  FREEI (*  ’ .NOLD. 1) 

CALL  FREEI ( ' I ' , LMOLD (1 ) , 2) 

IF (TYPE. EQ. ’RES’ ) THEN 
CALL  FREER  CR’  ,R.  1) 

CALL  FREER (’A* , A, 1) 

ELSEIF (TYPE. EQ. ’FLO' ) THEN 
CALL  FREER (’W’ ,R, 1) 

CALL  FREER (’C’ .A. 1) 

ENDIF 


CALL  RESIS0 (ID. NOLD. LMOLD . K . R . A , NNOD . NEQN . MBAN . MAXBAN. TYPE.  ERR) 
C 

C— 3.0  GET  NEXT  LINE  OF  DATA 
C 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END" 

CALL  FREEH (’  '. ENDFLAG. 3 . 1 ) 

IF  (ENDFLAG. EQ.  ’END’ ) THEN 
RETURN 
ENDIF 

C GET  NEW  ELEMENT  INFORMATION 

CALL  FREEI (*  ’.NNEW.l) 

CALL  FREEI  ('X'.LMJEWU)  ,2) 

CALL  FREEI ('N’ . INCR. 1) 

IF (INCR . EQ. 0)  INCR=1 
IF (TYPE. EQ. 'RES’ ) THEN 
CALL  FREER ('R' ,R. 1) 

CALL  FREER (’A' , A. 1) 

ELSEIF (TYPE. EQ. ’FLO’)  THEN 
CALL  FREER  <’W’ ,R, 1) 

CALL  FREER ('C' , A, 1) 

ENDIF 

C CHECK  NUMERICAL  ORDER 

IF (NNEW . LE . NOLD)  THEN 
WRITE (NTM. 2300)  NNEW 

WRITE (MOT. 2300)  NNEW 

ERR- . TRUE . 

RETURN 

ENDIF 

C OTNERATE  MISSING  ELEMENTS 

IF (NNEW . GT . NOLD+ 1 ) THEN 
DO  34  N—NOLD+ 1 , NNEW- 1 , 1 
DO  32  I-l.NDOF 

32  LMOLD (I)  - LMOLD (I)  ♦ INCR 

34  CALL  RESIS0 (ID. N. LMOLD. K.R. A, NNOD. NEQN. MBAN. MAXBAN, TYPE. ERR) 
ENDIF 

C DO  NEW  ELEMENT 

NOLD  - NNEW 
DO  36  I -1, NDOF 
36  LMOLD  (I)  - LWEW(I) 

CALL  RESIS0 (ID. NOLD. LMOLD. K. R. A. NNOD. NEQN. MBAN. MAXBAN. TYPE.  ERR) 


D 


4 


GO  TO  30 


NDOF  - 1 


2300  FORMAT ('  *•**  ERROR:  Element  number  ’.IS,*  Is  out  of  order.*) 
END 


RESISO 

SUBROUTINE 

RESISO (ID, NELM. LM.K. R, A. NNOD, NEQN . KBAN . MAXBAN . TYPE . ERR) 

C-SUB : RESISO  - REPORTS  ELEM  INFORMATION.  CHECKS  BANDWIDTH, 

C FORMS  ELEM  ARRAYS  4 ADDS  THEM  TO  SYS  ARRAYS 

C 

IMPLICIT  REAL *8  (A-H.O-Z) 

C 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR, NIN . NOT, ND1 , ND2 , ND3 , ND4 
C 

C RESISO:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 K (NEQN . MBAN ) . S (2, 2) 

INTEGER  ID (NNOD),  LM(2) 

LOGICAL* 1 ERR 
CHARACTER  TYPE* 3 
SAVE  /IOLIST/ 

C 

C REPORT  ELEMENT  INFORMATION  TO  OUTPUT  DATA  FILE 

C 

WRITE (NOT. 2000)  NELM.LM(l) ,LM(2) ,R.A 
C 

C NODE  ERROR  TRAP 

C 

DO  200  N=1 . 2 
NN*LM(N) 

IF (NN.LE. O.OR.NN.GT.NNOD)  THEN 
WRITE  (NTM.  2010)  NN 

WRITE  (NOT.  2010)  NN 
ERR= . TRUE . 


RETURN 

200  ENDIF 
C 

C DETERMINE  ELEMENT  BANDWIDTH 

C 

CALL  MBAND (ID, LM, 2. NNOD. MBAN, MAXBAN, ERR) 

C 

C FORM  ELEMENT  ARRAYS 

C 

IF (TYPE. EQ. ’FLO*)  THEN 
A = R*A 

R = 1.0 

ENDIF 

CALL  RESIS (S, R. A) 

C 

C ADD  ELEMENT  CONTRIBUTION  TO  SYSTEM  ARRAYS 

C 

CALL  ADDCM (ID. LM. S , K, 2 , 2 . NNOD, NEQN , MBAN) 


C 

C— 1.0  WRITE  ELEMENT  HEADER 
C 

WRITE (NOT,  2100) 

2100  FORMAT ( 

./*  CAPACITANCE  ELEMENTS’// 

Elere  I-Node  Mass  Capac.') 

C 

C — 2.0  GET  FIRST  ELEMENT.  FORM  4 ADD  ELEMENT  TO  SYS 
C 

INCR  « 0 

MASS  « 0.0 

CAP  - 0.0 

CALL  FREE 

CALL  FREETY 

CALLFREEK*  *.NOLD,l) 

CALL  FREEI ( * I ’ . LMOLD , 1 ) 

CALL  FREER ( *M’ .MASS, 1) 

CALL  FREER ( 'C* .CAP, 1) 

CALL  CAP AC 0 (ID. NOLD. LMOLD. C, MASS. CAP. NNOD. NEQN, MBAN, ERR) 


C 

C— 3.0  GET  NEXT  LINE  OF  DATA 
C 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH (’  ’ .ENDFLAG. 3. 1) 
IF (ENDFLAG.EQ. ’END’ ) THEN 
RETURN 
ENDIF 

C GET  NEW  ELEMENT  INFORMATION 

CALLFREEI(’  *,NNEW,1) 

CALL  FREEI ( ' I ' , LMNEW, 1 ) 

CALL  FREEI  CN’ . INCR.  1) 

IF (INCR . EQ. 0)  INCR-1 
CALL  FREER  CM’ .MASS,  1) 

CALL  FREER  CC' .CAP,  1) 

C CHECK  NUMERICAL  ORDER 

IF (NNEW.LE.NOLD)  THEN 
WRITE  (NTM. 2300)  NNEW 


WRITE (NOT, 2300)  NNEW 


ERR= . TRUE . 


RETURN 

ENDIF 

C GENERATE  MISSING  ELEMENTS 

IF (NNEW . GT . NOLD+1)  THEN 
DO  34  N*NOLD+l , NNEW-1 , 1 
LMOLD  « LMOLD  ♦ INCR 

34  CALL  CAPAC0 (ID,N. LMOLD. C. MASS. CAP. NNOD, NEQN. MBAN.  ERR) 
ENDIF 

C DO  NEW  ELEMENT 

NOLD  - NNEW 
LMOLD  * LMNEW 

CALL  CAP AC 0 ( ID. NOLD, LMOLD. C , MASS . CAP , NNOD , NEQN. MBAN . ERR) 


RETURN 


GO  TO  30 


2000  FORMAT (3 (5X, 15) .4X.2G10.3) 

2010  FORMAT (’  ••**  ERROR:  (Generated)  Node  *,15.*  Is  out  of  range.’) 
END 

C RESIS 

SUBROUTINE  RESIS (XK, R, A) 

C SUB: RESIS  - 2-NODE  CONVENTIONAL  RESISTANCE  ELEM 

C 

IMPLICIT  REAL *8  (A-H.O-Z) 

C 

DIMENSION  XK (2,2) 

C FORM  2X2  CONDUCTANCE  MATRIX 

XK  (1,1)  - A/R 
XK  (1,2)  - -A/R 
XK  (2,1)  - -A/R 
XK  (2,2)  - A/R 
RETURN 
END 

CAPACI 

SUBROUTINE  CAPACI (ID. C. NNOD. NEQN . MBAN. ERR) 

C-SUB: CAPAC I - FORMS  4 ASSEMBLES  CAPACITANCE  ELEMENTS 
C 

IMPLICIT  REAL*  6 (A-H.O-Z) 

C 

• c 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR. NIN . NOT. ND1 . ND2 . ND3 . ND4 
C 

C CAPACI:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 C (NEQN. MBAN) . MASS.  CAP 
INTEGER  ID (NNOD).  L8tJEW,  LMOLD 
LOG I CAL *1  ERR 
CHARACTER  ENDFLAG* 3 

SAVE  NNEW, NOLD, LHJEW, LMOLD. MASS . CAP . /IOLIST/ 


2300  FORMAT ( ’ ••••  ERROR:  Element  number  '.15. * Is  out  of  order.') 

END 

C -CAPAC0 

SUBROUTINE  CAPAC0 (ID . NELM. LM, C. MASS . CAP . NNOD. NEQN , MBAN. ERR) 

C-SUB :CAPAC0  - REPORTS  ELEM  INFORMATION. 

C FORMS  ELEM  ARRAYS  4 ADDS  THEM  TO  SYS  ARRAYS 

C 

IMPLICIT  REAL* 8 (A-H.O-Z) 

C 

C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR . NIN , NOT. ND1 , ND2 . ND3 . ND4 
C 

C FORM:  DATA  4 COMMON  STORAGE 

C 

REAL *8  C (NEQN, MBAN) , HASS,  CAP,  HC 
INTEGER  ID (NNOD).  LM  .NELM 
LOG I CAL *1  ERR 
SAVE  /IOLIST/ 

C 

C REPORT  ELEMENT  INFORMATION  TO  OUTPUT  DATA  FILE 

C 

WRITE (NOT. 2000)  NELM, LM, MASS, CAP 
C 

C NODE  ERROR  TRAP 

C 

NN-LM 

IF (NN.LE. O.OR.NN.GT.NNOD)  THEN 
WRITE  (NTM, 2010)  NN 

WRITE (NOT, 2010)  NN 
ERR*. TRUE. 

RETURN 

ENDIF 


D -5 


FORM  CLEMENT  ARRAY 


PR (2)  - 0.0D0 


C 

c 

c 

c 


c 

c 


HC  - MASS  • CAP 

ADD  ELX*»JT  CONTRIBUTION  TO  SYSTEM  ARRAYS 
CALL  ADOCM(ID.LM.MC.C.1,1.NNOD,NEQN.MBAN) 


RETURN 


PR  (3)  - 1.0D0 
ELSEIF  (NDOF.EQ.4)  THEN 
PR  (1 ) - -l.ODO 
PR  (2)  - -1.0DO/3.0DO 
PRO)  - 1.0D0/3.0D0 


2000  FORMAT  (2  (5X, 15) . 4X, 2G10.3) 

2010  FORMAT  ( * ••••  ERROR:  (Generated)  Node  '.15.'  Is  out  of  range.') 
END 


C 

C— 1.0  WRITE  ELEMENT  HEADER 
C 

WRITE  {KOT,  2100) 

2100  FORMAT ( 

./•  «=  ISOPARAMETRIC  ID  2 TO  4 -NODE  CONDUCTION  ELEMENTS'// 

Kx  « Conductivity  P * Density  C = Capacity'/ 

A * Area  Q = Int . Heat’// 

.*  Ele®  *©de  Coords.  Kx  P C A Q’ 

.*  Q»j*d . Points') 

C 

C — 2.0  GET  FIRST  ELEMENT,  FORM  £ ADD  ELEMENT  TO  SYS 
C 

I NCR  * 0 
DO  10  1*1. 1 
10  PX(I)  « 0.0 
PK  « 0.0 
PP  = 0.0 
PC  = 0.0 
PA  * 0.0 
PQ  » 0.0 
CALL  FREE 
CALL  F1JEETY 
CALL  FREE  I ( ' '.MOID.l) 

DO  15  1*1.4 
15  LMOLD (I)  « 0 

CALL  FREEI  ('I' ,LMOLD(l)  ,4) 

WDOr  * 4 

IF(LMOLD(4) .LZ.0)  N DOF-3 
SF(LMOLD(3) .LE.0)  WDOF-2 
CALL  FREER  ( * X ' , PX  ( 1 ) . 4 ) 

CALL  FREER  ( ' K * ,PK,  1) 

CALL  FREER  (*P* , PP,  1) 

CALL  FREER  ('C.  PC.  1) 

CALL  FREER  ( 'A*  , PA,  1) 

CALL  FREEACQ'.FQ.l) 

QUADFLAG  ■ *12345’ 

CALL  FREEH  CR\ QUADFLAG.  5.1) 

GAUSS  - .FALSE. 

IF  (QUADFLAG. EQ.  'GAUSS')  GAUSS  - .TRUE. 

IF  ( .WOT. GAUSS)  THEN 
C SET  DEFAULT  QUAD  POINTS 

IF  ‘WDOF.EQ.2)  TfiEN 


PR  ( 4 ) - l.ODO 
ENDIF 

C GET  QUAD  POINTS  FOR  EVALUATION  OF  [C ) 

CALL  FREER <’R' .PR (1) .NDOF) 

DO  20  1*1. NDOF 

IF  (ABS (PR(I) ) .GT. 1 .0D0)  THEN 

WRITE (NTH. 2200) 

WRITE  (NOT.  2200) 

20  ENDIF 
ENDIF 

2200  FORMAT ( ' WARNING:  Quadrature  points  should  normally 

♦ be  In  the  rar^ge  of  -1.0  to  **-1.0') 

CALL  ISO1D0 (ID. NOLD. LMOLD.C, K. EO. PX. PK. PA. PP. PC. PQ. PR, 
4NDOF . GAUSS , NNOO . NEON . KB AN . MAX BAN . ERR ) 

C 

C--3.0  GET  NEXT  LINE  OF  DATA 
C 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH  ( ' ' . ENDFLAG. 3.1) 

IF (ENDFLAG . EQ. ' END ' ) THEN 
RETURN 
ENDIF 

C GET  NEW  ELEMENT  INFORMATION 

CALLFREEIC  '.NNEW.l) 

DO  35  1=1.  4 
35  LMNEW(I)  = 0 

CALL  FREEI ( ' I ’ . LMNEW (1 ) , 4) 

NDOF  = 4 

IF (LMNEW (4) .LE. 0)  NDOF=3 
IF (LMNEW (3)  .LE.0)  NDOF =2 
CALL  FREEI ('N' . INCR. 1) 

IF(INCR.EQ.O)  INCR=1 
CALL  FREER ('X' . PX(1) , 4) 

CALL  FREER ('K' . PK.  1) 

CALL  FREER ('P' .PP.  1) 

CALL  FREER  CC'  .PC.  1) 

CALL  FREER ('A' ,PA,1) 

CALL  FREER ('Q' .PQ.  1) 

QUADFLAG  - '12345' 

CALL  FREEH ('R' .QUADFLAG. 5. 1) 

GAUSS  = .FALSE. 

IF (QUADFLAG. EQ. ’GAUSS')  GAUSS  = .TRUE. 

IF (.NOT. GAUSS)  THEN 
C SET  DEFAULT  QUAD  POINTS 

IF (NDOF. EQ. 2)  THEN 


PR  (1) 

* 

- SQRT (2 . 0D0/3 . 0D0) 

PR  (2) 

- 

- PR(1) 

ELSEIF  (NDOF. EQ. 3)  THEN 

PR  (1) 

- 

-1 . ODO 

PR  (2) 

- 

0.0D0 

PRO) 

- 

l.ODO 

ELSEIF (NDOF.EQ.4)  THEN 

PR(1) 

- 

-1 . ODO 

PR  (2) 

- 

-1.OD0/3.OD0 

PRO) 

- 

1 . 0D0/3 . ODO 

PR  (4) 

- 

l.ODO 

ENDIF 

C GET 

QUAD  POINTS  FOR  EVALUATION  OF  (C) 

CALL  FREER ('R' .PR (1) .NDOF) 


ISOID 

SUBROUTINE  ISOID (ID. XYZ.C, K. EO. NNOD, NEQN, MBAN. MAXBAN , ERR) 

C-SUB: ISOID  - FORMS  £ ASSEMBLES  ISOPARAM.1D  2-4-NCOE  FINITE  ELEMENTS 
C 

IMPLICIT  REAL* 8 (A-H.O-Z) 

C 

C 

C CAL -SAP:  DATA  £ COMMON  STORAGE 

C 

COMMOM  /IOLIST/  NTM . NTR . NIN . NOT . ND1 , ND2  , ND3  . ND4 
C 

C ISOID:  DATA  £ COMMON  STORAGE 

C 

C DICTIONARY  OF  VARIABLES  

C 

C VARIABLE 

C PK 

C PA 

C PP 

C PC 

C PQ 

C PX (4) 

C PR  (4) 

C*  ••**•***••••' 

REAL*  8 K (NEQN. MBAN)  . C (NEQN , MBAN ) ,XYZ (NNOD, 3) , EO (NNOD) . 

4 PK. PA. PP . PC, PQ,PX(4) ,PR(4) 

INTEGER  ID  (NNOD).  L«4EW(4),  LMOLD  ( 4 ) 

LOGICAL*  1 ERR.  GAUSS 
CHARACTER  ENDFLAG* 3.  QUADFLAG* 5 

SAVE  RSEW. NOLD. LWiEW, LMOLD, PK, PA.PP. PC. FQ. PX, PR, GAUSS. /IOLIST/ 


DESCRIPTION 

ELEMENT  CONDUCTIVITY 
AREA  AVAILABLE  FOR  HEAT  TRANSFER 
ELEMENT  DENSITY 
ELEMENT  CAPACITY 

ELEMENT  INTERNAL  HEAT  GENERATION  RATE/LENGTH 
ELEMENT  NODE  LOCAL  COORDINATES 
QUAD  POINTS  -1.0  TO  1.0 


PR  (1)  «=  - SORT  (2. 0D0./3. 0D0) 


DO  40  1=1. NDOF 


PR  (2)  * - PRO  ) 


IF  (ABS  (PR  (I) ) .GT. l.ODO)  THEN 


EL-SEIF  (NDOF . EQ.  3)  THEN 


WRITE (NTM,  2200) 


PRO)  * -1-000 


WRITE (NOT.  2200) 
40  ENDIF 


D 


6 


ENDIF 


C 

IS0241 


C CHECK  NUMERICAL  ORDER 

IF  (NNEW. LE. NOLD)  THEN 
IfRITE  (NTH,  2300)  NNEW 

WRITE  (NOT.  2300)  NNEW 


ERR= . TRUE . 


RETURN 

ENDIF 

C GENERATE  MISSING  ELEMENTS 

IF (NNEW .GT . NOLD+1 ) THEN 
DO  50  N-NOLD+l.NNEW-1,1 
DO  45  I-l.NDOF 

45  LMOLD(I)  - LMOLD (I)  ♦ INCR 

50  CALL  I SOI DO (ID. N , LMOLD, C . K . EO, PX, PK. PA, PP , PC, PQ, PR, 

+ NDOF. GAUSS. NNOD. NEON. MBAN. MAXBAN. ERR) 

ENDIF 

C DO  NEW  ELEMENT 

NOLD  - NNEW 
DO  55  I-l.NDOF 
55  LMOLD (I ) - LMNEW(I) 

CALL  ISO1D0 (ID. NOLD , LMOLD. C,  K, EO. PX . PK, PA, PP , PC. PQ. PR, 
♦NDOF,  GAUSS,  NNOD.NEQN.  MBAN,  MAXBAN  , ERR) 


SUBROUTINE  IS02  41  (KE. CE. EE.  PK, PA. PP. PC. PQ. PX , NELNOO, RPT, GAUSS ) 
C--SUB : IS0241  - 2 TO  4 NODE  ISOPARAMETRIC  CONDUCTION  FINITE  ELEMENT 


NODE  NUMBERING 


C DICTIONARY  OF  VARIABLES 


C VARIABLE 

C DUMMY  VARIABLES 
C NELNOO 

C NIP 

C KE (NELNOD.  NELNOD) 

C CE (NELNOD. NELNOD) 

C EE  (NELNOD) 

PK 


PA 

PP 

PC 

PQ 

PX  (NELNOD) 
RPT  ( 4 ) 


DESCRIPTION 

NUMBER  OF  ELEMENT  NODES  (2  TO  4] 

NUMBER  OF  INTEGRATION  POINTS 

ELEMENT  CONDUCTANCE  MATRIX 

ELEMENT  CAPACITANCE  MATRIX 

ELEMENT  INTERNAL  GENERATION  RATE  VECTOR 

CONDUCTIVITY 

AREA  AVAILABLE  FOR  HEAT  TRANSFER 
WEIGHT  DENSITY 
HEAT  CAPACITY 

INTERNAL  GENERATION  RATE  PER  UNIT  VOL. 
(GLOBAL)  COORDINATES  OF  ELEM.  NODAL  POINTS 


QUAD.  POINTS  USED  TO  EVALUATE  [C] 
C GAUSS 


GO  TO  30 

2300  FORMAT ('  BRROR : Element  number 


QUADRATURE  METHOD  TO  EVALUATE  (C) 

Is  out  of  order.')  C 


END 

C ISO1D0  C 

SUBROUTINE  ISO1D0 (ID . NELM. LM. C, K. EO. PX, PK. PA.PP. PC. PQ.  PR. 

♦NDOF. GAUSS. NNOD.NEQN. MBAN. MAXBAN. ERR) 

C-SUB : ISO1D0  - REPORTS  ELEM  INFORMATION.  CHECKS  BANDWIDTH. 

C FORMS  ELEM  ARRAYS  4 ADDS  THEM  TO  SYS  ARRAYS 

C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

c 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR, NIN , NOT, ND1 , ND2 . ND3 , ND4 
C 

C ISO1D0 : DATA  4 COMMON  STORAGE 

C 

REAL* 8 C (NEON. MBAN) . K (NEQN . MBAN ) .EO(NNOD) . CE (4 . 4 ) . KE (4 . 4 ) . 

♦EE  (4) . PX ( 4 ) .PR (4) 

INTEGER  ID (NNOD).  LM(4) 

LOGICAL* 1 ERR. GAUSS 
SAVE  /IOLIST/ 

C 

C REPORT  ELEMENT  INFORMATION  TO  OUTPUT  DATA  FILE 

C 

IF (GAUSS)  THEN 

WRITE  (NOT. 2000)  NELM, LM (1 ) ,PX(1) , PK, PP . PC. PA. PQ 

WRITE  (NOT.  2005)  (LM(I).PX(I)  ,1-2.  NDOF) 

ELSE 

WRITE (NOT. 2010)  NELM, LM (1 ) ,PX(1) . PK. PP . PC , PA. PQ. 

♦ (PR(I) . I-l.NDOF) 

WRITE  (NOT. 2005)  (LM (I ) . PX  (I ) . 1=2 . NDOF) 

ENDIF 

2000  FORMAT (215. 6G8.2E1. ’GAUSS’) 

2005  FORMAT (5X.I5.G6.2) 

2010  FORMAT (215, 6G8 . 2E1 . 4F6. 3) 

C 

C NODE  ERROR  TRAP 

C 

DO  200  N-l.NDOF 
NN-LM (N) 

IF  (NN.LE. O.OR.NN.GT.NNOD)  THEN 
WRITE  (NTM.  2020)  NN 

WRITE  (NOT,  2020)  NN 
ERR-. TRUE. 

RETURN 

200  ENDIF 

2020  FORMAT ('  ERROR:  (Generated)  Node  '.IS,'  Is  out  of  range.') 

C 

C DETERMINE  ELEMENT  BANDWIDTH 

C 

CALL  MBAND  (ID.  LM.NDOF,  NNOD,  MBAN.  MAXBAN,  ERR) 

C 

C FORM  ELEMENT  ARRAYS 

C 

CALL  IS0241 (KE, CE. EE, PK, PA.PP. PC, PQ.PX, NDOF, PR, GAUSS) 

C 

C ADD  ELEMENT  CONTRIBUTION  TO  SYSTEM  ARRAYS 

C 

CALL  ADDCM (ID, LM, KE, K, NDOF. 4 . NNOD, NEQN . MBAN ) 

CALL  ADDCM ( ID, LM.CE.C, NDOF, 4. NNOD. NECK. MBAN) 

DO  300  N-l.NDOF 
300  EO  (LM  (N) ) = EE  (N) 


RETURN 

END 


- .TRUE.  USE  GAUSS  QUADRATURE 


- .FALSE.  USE  GIVEN  QUAD  POINTS 
C INTERNAL  VARIABLES 
C QPT  ( 4 ) 

GAUSS  QUAD.  POINTS 
C QWT  (4) 

GAUSS  QUAD.  WEIGHTS 
C H ( 4 ) 

ELEMENT  SHAPE  FUNCTION  VALUES 
C DEL  (4) 

ELEMENT  LOCAL  DERIVATIVE  VALUES 
C R LOCAL  COORDINATE 

C JACOB  JACOBIAN 

C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

REAL*  8 KE (4. 4) . CE (4. 4) . EE (4) . PX (4)  , RPT(4) 

REAL*8  JACOB,  QPT(4),  QWT(4),  H (4) . DEL  ( 4 ) 

LOGICAL* 1 GAUSS 

C-1.0  ZERO  ELEMENT  ARRAYS 

DO  10  1=1.4 
EE  (I ) = 0.0D0 
EX)  10  J=1 , 4 
KE(I.J)  = 0.0D0 
CE(I.J)  = 0.0D0 
10  CONTINUE 

C-2.0  EVALUATE  CONDUCTANCE  MATRIX  AND  INTERNAL  GENERATION  RATE  VECTOR 
C-2.1  GET  GAUSS  QUAD  POINTS  AND  WEIGHTS 

C NOTE:  USE  NIP-NELNOD-1  FOR  EXACT  INTEGRATION  OF  [KE]  4 (EE) 

NIP  - NELNOD  - 1 

CALL  GAUSCO (NIP , QPT , QWT) 

DO  30  IP-1, NIP 

C-2.2  GET  SHAPE  FUNCTION  VALUES  AT  QUAD  POINT 
CALL  SHP241 (QPT (IP) ,H. NELNOO) 

C-2.3  GET  LOCAL  DERIVATIVE  VALUES  AT  QUAD  POINT 
CALL  DER241 (QPT (IP) .DEL, NELNOD) 

C-2.4  COMPUTE  JACOBIAN  AT  THE  POINT  (NOTE:  FOR  ID  JACOBIAN  IS  1X1) 
CALL  JACOB 1 (NELNOD . DEL . PX . JACOB ) 

C-2.5  ACCUMULATE  QUADRATURE  CONTRIBUTION  TO  ELEMENT  MATRICES 
DO  20  1-1 . NELNOO 

C ELEMENT  INTERNAL  GENERATION  RATE  VECTOR 

EE (I)  - EE (I)  ♦ QWT (IP) *PQ*H (I) ‘JACOB 
DO  20  J-l, NELNOO 

C ELEMENT  CONDUCTANCE  MATRIX 

KE(I.J)  - KE(I.J)  ♦ QWT  (IP)  *DEL  (I)  *PK*DEL  ( J)  / JACOB 
20  CONTINUE 
30  CONTINUE 

C-3.0  EVALUATE  CAPACITANCE  MATRIX 
C-3.1  GET  QUAD  POINTS  AND  WEIGHTS 

C NOTE:  USE  NIP^ELNOD  FOR  EXACT  INTEGRATION  OF  (CJ 
NIP  - NELNOO 

C USE  GAUSS  QUADRATURE  IF  SPECIFIED 

C EVALUATE  HEIGHT  COEFFICIENTS  IF  QUAD  POINTS  ARE  SPECIFIED 
ir (GAUSS)  CALL  GAUSCO (NIP, RPT, QWT) 

IF (.NOT. GAUSS)  CALL  QUADCO (NIP . QWT . RPT) 

DO  50  IP-1, NIP 

C-3.2  GET  SHAPE  FUNCTION  VALUES  AT  QUAD  POINT 
CALL  SHP241 (RPT (IP) ,H, NELNOO) 

C-3.3  GET  LOCAL  DERIVATIVE  VALUES  AT  QUAD  POINT 
CALL  DER241  (RPT (IP) .DEL, NELNOD) 

C-3.4  COMPUTE  JACOB I AN  AT  THE  POINT  (NOTE:  FOR  ID  JACOBIAN  IS  1X1) 
CALL  JACOBI  (NELNOD. DEL. PX. JACOB) 

C-3.5  ACCUMULATE  QUADRATURE  CONTRIBUTION  TO  ELEMENT  MATRICES 
DO  40  I -1, NELNOD 
DO  40  J-l, NELNOD 

CE  (I,  J)  - CE(I.J)  ♦ QWT(IP)  *H  (I)*PP*PC*H(J) ‘JACOB 


D -7 


40  CONTINUE 
50  CONTINUE 


RETURN 

END 


C 

SHP241 


SUBROUTINE  SHP241  (R.H, NELNOD) 

C-SUB : SRP24L  - DETERMINES  THE  VALUE  OF  2-4  NODE  ISOPARAMETRIC  ID 
C SHAPE  FUNCTIONS,  H,  AT  LOCAL  COORDINATE  R 


C 

C 


NODE  NUMBERING  1 3 4 2 

DICTIONARY  OF  VARIABLES 


C VARIABLE  DESCRIPTION 

C R LOCAL  COORDINATE  -1  TO  +1 

C H(2  TO  4)  VALUE  OF  SHAPE  FUNCTIONS 

C NELNOD  NUMBER  OF  ELEMENT  NODES  l 2 TO  4 ] 


IMPLICIT  REAL* 8 (A-H.O-Z) 
DIMENSION  H (NELNOD) 


H (1 ) - 0.5*(1.0  - R) 

3(2)  - 0.5*  (1.0  ♦ R) 

IF (NELNOD. GE. 3)  THEN 

B (1 ) » H (1)  - 0 . 5*  (1 . 0 - R*R) 


H<2)  ~ H (2)  - 0.5M1.0  - R*R) 


fl(3>  » (1.0  - R*R) 

ENDIF 

IF (NELNOD. EQ. 4)  THEN 

3(1)  « H (1)  ♦ (-9 . 0*R*R*R  R*R  ♦ 9.0*R  - 1.0)/16.0 


3(2)  « 3(2)  > (9. 0*R*R*R  + R*R  - 9.0*R  - 1.0)/16.0 

H (3 ) - H (3)  + (27 . 0*R*R*R  + 7.0*R*R  - 27.0*R  - 7.0)/16.0 

3(4)  « (-27 . 0 *R*R*R  - 9.0*R*R  ♦ 27.0*R  + 9.0J/16.0 
ENDIF 
RETURN 
END 


C 

DER241 


C-SUB 

C 

c 

c 

c 

c~ 

c 


SUBROUTINE  DER2 4 1 (R. DEL, NELNOD) 

:DER241  - DETERMINES  THE  VALUE  OF  LOCAL  DERIVATIVES.  DEL.  OF 
2-4  NODE  ISOPARAMETRIC  ID  SHAPE  FUNCTIONS  AT  LOCAL  COORDINATE  R 

NODE  NUMBERING  1—3 4 2 

DICTIONARY  OF  VARIABLES  

VARIABLE  DESCRIPTION 

R LOCAL  COORDINATE  -1  TO  +1 

DEL (2  TO  4)  VALUE  OF  LOCAL  DERIVATIVE  OF  SHAPE  FUNCTIONS 
NELNOD  NUMBER  OF  ELEMENT  NODES  [ 2 TO  4 ) 


IMPLICIT  REAL*  8 (A-H.O-Z) 
DIMENSION  DEL (NELNOD) 


DEL (1 ) - -0.5 
DEL (2)  * 0.5 
IF (NELNOD. GE. 3)  THEN 
DEL (1 ) - DEL (1)  + R 


REAL*0  JACOB 

DIMENSION  DEL  (NELNOD).  XCOORD (NELNOD) 
JACOB  - 0.0 
DO  10  N«l. NELNOD 

10  JACOB  - JACOB  ♦ DEL (N) • XCOORD (N) 
RETURN 
END 


C 

GAUSCO 


SUBROUTINE  GAUSCO (N , A. W) 

C— SUBrGAUSCO  - GAUSS  QUADRATURE  ABSCISSAE  AND  WEIGHT  COEFFICIENTS 
C 

C DICTIONARY  OF  VARIABLES  


VARIABLE 

N 

A (N) 

M(N) 


DESCRIPTION 

NUMBER  OF  QUADRATURE  POINTS 

ABSCISSAE 

WEIGHTS 


IMPLICIT  REAL*0 (A-B.O-Z) 
DIMENSION  A (N) . W(N) 
IF(N.EQ.l)  THEN 
A (1)  » 0.0D0 


W(l) 

- 

2.0DC 

ELSEIF (N .EQ. 2)  THEN 

A (1)  = -1.0D0/SQRT(3. 0D0) 

A (2) 

- 

"A  (1) 

W(l) 

- 

1.0D0 

W (2) 

1 . 0D0 

ELSEIF (N.EQ. 3)  THEN 
A (I)  - -SQRT (3. 0/5.0) 

A (2) 

- 

0.OD0 

A (3) 

- 

-A(l) 

W(l) 

- 

5. 0D0/9.0D0 

W(2) 

- 

8. 0D0/9.0D0 

N(3) 

. 

W(l) 

ELSEIF (N.EQ. 4)  THEN 

A (1)  - -0.861136311S94053D0 

A(2) 

« 

-0. 33998104358485 6D0 

A (3) 

= 

-A  (2) 

A (4) 

= 

-A(l) 

W(l) 

- 

0.347854845137454D0 

W(2) 

- 

0 . 6521 451 54 862546D0 

W(3) 

= 

W(2) 

W (4) 

B 

W(l) 

ELSE 

PAUSE  ' ***•  ERROR: SUB: GAUSCO:  Gauss  coefs.  not  available.’ 
ENDIF 
RETURN 
END 


DEL (2)  « DEL  (2)  ♦ R 

DEI.  (3)  « -2  0*R 

ENDIF 

IF  (NELNOD. EQ. 4)  THEN 

DEL (1)  - DEL (1)  + (-27 . 0*R*R  4 2.0*R  * 9.0)/16.0 

DEL (2)  *=  DEL  (2)  ♦ (27.0*R*R  + 2.0*R  - 9.0)/16.C 

D£L(3)  * DEL (3)  + <81.0*R*R  + 14.0*R  - 27.0) /16.0 

DEL (4)  « (-81 *R*R  - 18*R  ♦ 27.0)/16.0 
ENDIF 
RETURN 
END 


C 

QUADCO 


SUBROUTINE  QUADCO (N.W.R) 

C--SUB :QUADCO  - EVALUATES  WEIGHT  COEFS.  GIVEN  QUAD  POINTS 
C 

C DICTIONARY  OF  VARIABLES  

C 

C VARIABLE  DESCRIPTION 

C INPUT 

C N NUMBER  OF  QUADRATURE  POINTS 

C R (N)  SPECIFIED  QUAD  POINTS 

C OUTPUT 

C W(N)  WEIGHT  COEFFICIENTS 

C LOCAL 
C C (4 , 4) 


JACCBl 

SUBROUTINE  JACOBI (NELNOD . DEL . XCOORD , JACOB) 
e-SUS: JACOBI  - COMPUTES  THE  JACOBIAN,  JACOB,  FOR  2-4  NODE  ISOPARA. 

C E LJSK.  USING  VALUES  OF  LOCAL  DERIVATIVES  OF  SHAPE  FUNCTION  IN  DEL 

C 

6 » J « dx/dr  - d/dr  {HI  H2  H3  H4)(X1  X2  X3  X4)T  ; dEi/dr  * DELI 

e 

e KOBE  NUMBERING  1 3 4 2 

C DICTIONARY  OF  VARIABLES  

C 

C VARIABLE  DESCRIPTION 

C*  NELNOD  NUMBER  OF  ELEMENT  NODES  [2  TO  4) 

C DEL  (2  TO  4)  VALUE  OF  LOCAL  DERIVATIVE  OF  SHAPE  FUNCTIONS 

C XCOORD  (2  TO  4)  (GLOBAL)  COORDINATES  OF  NODES 

C JAC 

C 

ILLICIT  REAL*  8 (A-a.O-Z) 


QUADRATURE  RULE  EQUATION  COEFFICIENTS 
C 

IMPLICIT  REAL* 8 (A-H.O-Z) 

DIMENSION  W( 4) . R(4),  C(4,4> 

DO  10  I-l.N 

W(I)  - (1.0D0  - (-1 . 0D0) **I) /DBLE (I) 

DO  10  J*1,N 

10  C(I,  J)  - R ( J)  **  (1-1) 

CALL  CROUT (C,W,N,4,1,0) 

RETURN 

END 

C IS02D4 

SUBROUTINE  IS02D4 (ID. XYZ, C, K. EO, NNOD , NEQN , MB AN , MAXBAN , ERR) 

C-SUB :IS02D4  - FORMS  AND  ASSEMBLES  ISOPARAMETRIC  2D  4-NODE  ELEMENTS 
C 

IMPLICIT  REAL*  0 (A-H.O-Z) 


-8 


D 


IMPLICIT  REAL*  8 (A-H.O-Z) 


C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR. NIN , NOT. ND1 . ND2 . ND3 , ND4 
C 

C RESIST:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 K (NEQN, MBAN) . C (NEQN . MEAN ) .XYZ  CNNOD, 3) .EO(NNOD) ,PK(3) 
INTEGER  ID  (NNOD).  U«EW(4),  LMOLD (4) 

LOGICAL*  1 ERR 

CHARACTER  END FLAG*  3 , REG* 3 

SAVE  LMNEW, LMOLD, NNEW, NOLD. PK, PC, PP , PQ, PD, REG, /IOLIST/ 

NDOF  * 4 
C 

C— 1.0  WRITE  ELEMENT  HEADER 
C 

WRITE (NOT, 2100) 

2100  FORMAT ( 

./'  «*  ISOPARAMETRIC  2D  4 -NODE  CONDUCTION  ELEMENTS'// 

. * c - Capacity  P * Density  0 * Internal  Heat  Rate’/ 

. * D - Thickness  for  REG*PLN;  Subtended  angle  for  REG* AX I ’ / / 

. * Elem  I J K L Rn  Kyy  Kxy  C 

POD  REG’) 

C 

C — 2 . 0 GET  FIRST  ELEMENT.  FORM  4 ADD  ELEMENT  TO  SYS 
C 

PQ  • 0.0 
PD  * 1.0 
REG  * * PLN ’ 

INCR  « 1 
CALL  FREE 
CALL  FREETY 
CALLFREEIC  ’,NOLD,l) 

CALL  FREEI ( ’ I ’ . LMOLD ( 1 ) , 4 ) 

CALL  FREER ( ’ K ’ , PK (1 ) , 3 ) 

CALL  FREER ('C', PC. 1) 

CALL  FREER ( 'P ’ .PP.l) 

CALL  FREER  CQ’.  PQ.l) 

CALL  FREER  CD’,  PD,  1) 

CALL  FREEH ('G'. REG, 3,1) 

CALL  ISO2D0 (ID. XYZ . NOLD, LMOLD. K, C, EO. PK. PC. PP,  PQ. PD.  REG. 

. NNOD, NEQN . MBAN , MAX BAN , ERR ) 

C 

C — 3 . 0 GET  NEXT  LINE  OF  DATA 
C 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH ( ' ’ .ENDFLAG. 3. 1) 

IF(ENDFLAG.EQ. ’END’ ) THEN 
RETURN 
ENDIF 

C GET  NEW  ELEMENT  INFORMATION 

CALL  FREEI (’  ’ .NNEW.l) 

CALL  FREEI (’I’ .LMNEW(l) , 4) 

CALL  FREEI (’N’ . INCR, 1) 

IF (INCR .EQ. 0)  INCR=1 
CALL  FREER ('K' ,PK (1)  . 3) 

CALL  FREER (*C’ .PC. 1) 

CALL  FREER (’P’ .PP.  1) 

CALL  FREER  CQ’ .PQ.  1) 

CALL  FREER  CD’,  PD.  1) 

CALL  FREEH (’G’ .REG. 3, 1) 

C CHECK  NUMERICAL  ORDER 

IF  (NNEW. LE. NOLD)  THEN 
WRITE  (NTM,  2300)  NNEW 

WRITE  (NOT.  2300)  NNEW 


ERR- . TRUE . 


RETURN 

ENDIF 

C GENERATE  MISSING  ELEMENTS 

IF  (NNEW .GT . NOLD+1 ) THEN 
DO  34  N-NOLD+1. NNEW- 1.1 
DO  32  I-l.NDOF 

32  LMOLD (I)  * LMOLD (I)  ♦ INCR 

34  CALL  ISO2D0 (ID, XYZ , N. LMOLD, K, C, EO, PK, PC. PP, PQ. PD, REG. 

. NNOD . NEQN , MBAN . MAXBAN . ERR ) 

ENDIF 

C DO  NEW  ELEMENT 

NOLD  - NNEW 
DO  36  1*1, NDOF 
36  LMOLD (I)  * LMNEW(I) 

CALL  ISO2D0 (ID. XYZ . NOLD. LMOLD, K, C, EO. PK. PC, PP, PQ, PD. REG, 

. NNOD , NEQN , MBAN , MAXBAN , ERR ) 

GO  TO  30 

2300  FORMAT (’  ****  ERROR:  Element  number  ’.15.’  Is  out  of  order.') 
END 


c ISO2D0 

SUBROUTINE  ISO2D0 (ID. XYZ . NELM, LM.K. C. EO. PK. PC, PP. PQ, PD. REG. 

. NNOD . NEQN . MBAN . MAXBAN , ERR ) 

C-SUB : ISO2D0  - REPORTS  ELEM  INFORMATION.  CHECKS  BANDWIDTH, 

C FORMS  ELEM  ARRAYS  4 ADDS  THEM  TO  SYS  ARRAYS 

C 


c 

c 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR, NIN. NOT. ND1 . ND2 , ND3 . ND4 
C 

C FORM:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 XY Z (NNOO , 3 ) , K (NEQN . MBAN ) , C (NEQN. MBAN)  , EO (NNOD) . 
.PK<3) ,KE (4. 4) ,CE (4. 4) .EE (4) . XY (2. 4) 

INTEGER  ID (NNOD).  LM(4) 

LOGICAL* 1 ERR 
CHARACTER  REG* 3 
SAVE /IOLIST/ 

C 

C REPORT  ELEMENT  INFORMATION  TO  OUTPUT  DATA  FILE 

C 

WRITE (NOT. 2000)  NELM. LM. PK. PC, PP. PQ. PD. REG 
2000  FORMAT (6X.  14,  IX, 413, IX, 7G8 . 2, A3) 

C 

C NODE  ERROR  TRAP 

C 

DO  200  N-1,4 
NN-LM(N) 

IF (NN . LE . 0 . OR . NN . GT . NNOD ) THEN 
WRITE  (NTM, 2010)  NN 

WRITE (NOT. 2010)  NN 
ERR* . TRUE . 


RETURN 

200  ENDIF 
C 

C DETERMINE  ELEMENT  BANDWIDTH 

C 

CALL  MBAND (ID. LM. 4, NNOO. MBAN. MAXBAN. ERR) 

C 

C FORM  ELEMENT  COORD  SUB-ARRAY  XY ( 2 . 4 ) 

C 

DO  300  1*1,2 
DO  300  J*1 , 4 

300  XY(I.J)  « XYZ  (LM  ( J)  , I) 

C 

C FORM  ELEMENT  ARRAYS 

C 

CALL  IS02D (KE.CE, EE. PK. PC.PP, PQ. PD. REG. XY) 
C 

C ADD  ELEMENT  CONTRIBUTION  TO  SYSTEM  ARRAYS 

C 

CALL  ADDCM (ID. LM. KE, K. 4 , 4 .NNOD, NEQN. MBAN) 
CALL  ADDCM (ID, LM, CE, C, 4 . 4 .NNOD. NEQN . MBAN) 

DO  400  N=1 , 4 
400  EO  (LM  (N)  ) - EE  (N) 


RETURN 

2010  FORMAT ( ’ **•*  ERROR:  (Generated)  Node  *.15,’  Is  out  of  range.') 

END 

c IS02D 

SUBROUTINE  IS02D (KE, CE. EE.PK. PC. PP. PQ, PD. REG. XY) 

C-SUB :IS02D  - 2D  ISOPARMETRIC  4-NODE  CONDUCTION  FINITE  ELEMENT 
C MODIFICATION  OF  SUBROUTINE  DEVELOPED  BY  R.  TAYLOR 

C U.C.  BERKELEY  FOR  PROGRAM  *HEAT* 

C 

C DICTIONARY  OF  VARIABLES 

C 

C VARIABLE 

C KE  ( 4 , 4 ) 

C CE (4,4) 

C EE  ( 4 ) 

C PK  (3) 

C PC 

C PP 

C PQ 

C PD 

C REG 

C XY (2,4) 

C 

IMPLICIT  REAL* 8 (A-H.O-Z) 

CHARACTER  REG* 3 

REAL* 8 SG (4) ,TG(4), KE (4 , 4) , CE (4 , 4) , EE (4) , XY (2,4) , PK (3) 

COMMON  / ISODAT/SH (3 , 4) , DV 
SAVE  /ISODAT/ 

DATA  SG/1 . , 1 . . -1 . . -1 . /. TG/-1 . . 1 . , 1 . , -1 . / 

G « 1. 0/SQRT (3.0) 

DO  100  1-1,4 
EE (I)  - 0.0 
DO  100  J-1.4 
CE(I.J)  - 0.0 

100  KE(I.J)  * 0.0 
PCPP*PC*PP 
DO  103  L— 1 , 4 

CALL  I SOS HP (SG (L) *G, TG (L) *G. XY) 

RR  * 0.0 
DV  * DV*PD 
DO  101  1-1.4 

101  RR  * RR  + XY  (1 . 1) *SH (3,1) 

IF (REG.EQ. 'AXI ’ ) DV  * DV*RR 
DO  102  J-1,4 
SHJ  - SH(3,J)*DV 


DESCRIPTION 

ELEMENT  CONDUCTANCE  MATRIX 

ELEMENT  CAPACITANCE  MATRIX 

ELEMENT  INTERNAL  GENERATED  FLUX  VECTOR 

ELEMENT  CONDUCTIVITY  TENSOR:  Kxx,  Kyy.  Kxy 

ELEMENT  SPECIFIC  HEAT  CAPACITY 

ELEMENT  DENSITY 

VOLUMETRIC  INTERNAL  HEAT  GENERATION  RATE 
ELEMENT  THICKNESS /SUBTENDED  ANGLE 
ELEMENT  TYPE:  'PLN'  OR  ’AXI’ 

ELEMENT  NOOAL  COORDINATES 
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C TREE:  BCOND  - TREE  OF  SUBROUTINES  TO  ESTABLISH  BOUNDARY  CONDITIONS 


EE  ( J)  - EE(J)  ♦ SH J*PQ 
SHJ  - SRJ*PCPP 

A1  - (PK(1) *SH(1.  J)  ♦ PK  (3)  *SH  (2,  J)  ) *DV 
A2  - (PK  (3)  *SH  (1,  J)  ♦ PK  (2)  *SH  (2 , J)  ) *DV 
DO  102  1-1,4 

KE(I.J)  - KE(I.J)  ♦ A1  * SH  (1,1)  «■  A2*SH  (2,1) 

102  CE(I.J)  - CE(I.J)  ♦ SHJ*SH  (3,1) 

103  CONTINUE 
RETURN 
END 

C ISOSHP 

SUBROUTINE  ISOSHP (SS . TT, XY ) 

C-SUB: ISOSHP  - SHAPE  FUNCTION  SUBROUTINE  FOR  IS02D 
IMPLICIT  REAL*  8 (A-H.O-Z) 

REAL* 8 SX (2,2) ,XS (2. 2) ,R  (4) ,T(4) ,XY (2. 4) 

COMMON  /ISOOAT/SH (3.4) ,DV 
SAVE  /ISOOAT/ 

DATA  R/-0.5, 0.5. 0.5. -0 . 5/. T/-0 . 5, -0 . 5.  0 . 5.  0 . 5/ 

DO  100  1-1.4 

SH  (3, 1 ) - (0.5+R(I)*SS)* (0. 5+T (I) *TT) 

SHU, I)  - R (I ) • (0 . 5+T  (I) *TT) 

100  SH  (2, 1)  - T (I) • (0.5+R(I)*SS) 

DO  101  1-1,2 

DO  101  J-1.2 
XS  (I,  J)«0.0 
DO  101  K-1,4 

101  XS(Z.J)  - XS(I.J)  ♦ XT  (I.K)  *SH  (J.K) 

DV  « XS  (1, 1) *XS (2. 2)  - XS  (1 , 2) *XS (2.1) 

SX  (1,1)  - XS (2,2) /DV 
SX  (2. 2)  « XS (1.1) /DV 
SX  (1, 2)  - -XS (1,2) /DV 
SX  (2. 1)  - -XS (2,1) /DV 
DO  102  K-l.  4 

CC  - SH  (1 , K) * SX  (1,1)  + SH (2 , K) *SX (2 , 1) 

SH  (2,  K)  * SH  (1 , K)  *SX  (1,2)  ♦ SH  (2 . K)  *SX  (2 . 2) 

102  SH  ( 1 , K)  - CC 
RETURN 

END 

C VNMRT 

SUBROUTINE  VNMRT (ID, K, NN CO, NEON . MB AN. MAX BAN, ERR) 

C-SUB: VNMRT  - FORMS  AND  ASSEMBLES  VARIABLE  NODE  MRT  ELEMENTS 
C ••*  PRESENTLY  A STUB 
RETURN 
END 


C MB  AND 

SUBROUTINE  MBAND (ID, LM , NDOF , NNOD . MB AN , MAXBAN , ERR) 

C-SUB:  MBAND  - DETERMINES  ELEMENT  BAND  WIDTH 
C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 


c 

C RBCOND 

SUBROUTINE  RBCOND (MP ID, NNOO, ERR) 

C-SUB: RBCOND  - ROOT  SUBROUTINE  TO  BCOND 

C BCOND  - ESTABLISHES  BOUNDARY  CONDITION  FLAGS  BY  MODIFYING 

ID  (NNOD) 

C 


C DICTIONARY  OF  VARIABLES 

C 

C VARIABLE  DESCRIPTION 

C VARIABLES  PASSED  TO  SUBROUTINE 
C ERR  DO-WHILE  TERMINATOR  FLAG 

C NNOD  NUMBER  OF  NODES  IN  SYSTEM 

C MBAN  (HALF)  BANDWIDTH  OF  SYSTEM  EQUATIONS 

C MPID  POINTER  TO  ID (NNOD)  IN  BLANK  CO***ON 

C ID (NNOO)  EQUATION  NUMBER/B.C.  FLAG  ARRAY; 

C ABS (ID (N) ) - EQUATION  NUMBER  OF  NODE  N 

C SIGN (ID  (N) ) - -1  - TEMPERATURE  PRESCRIBED  NOOE 

C SIGN (ID  (N))  - +1  - FLUX  PRESCRIBED  NODE 


C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  MTOT.NP. IA(1) 

COMMON  /DBSYS/  NUMA,  NEXT. IDIR. IP (3) 

COMMON  /IOLIST/  NTM, NTR, NIN , NOT . ND1 , ND2 . ND3 , ND4 
C 

C RBCOND:  DATA  4 COHION  STORAGE 

C 

LOGICAL* 1 ERR 

SAVE  /DBSYS/, /IOLIST/ 

C 

C--0.0  WRITE  HEADER 
C 

WRITE (NTM, 2000) 

WRITE (NOT, 2000) 

2000  FORMAT (/'  «==  BOUNDARY  CONDITIONS  AND  EQUATION  NUMBERS’) 


C 

C--1.0  FIND  SEPARATOR  "BOUNDARY" 

C 

CALL  FINDN (’BOUNDARY  \ 8. KEY) 
IF (KEY . EQ. 1 ) THEN 
WRITE  (NTM, 2100) 

WRITE (NOT. 2100) 

RETURN 

ENDIF 

CALL  FREETY 


C 

c 

c 

c 


COMMON  /IOLIST/  NTM, NTR, NIN, NOT, ND1 , ND2 , ND3 , ND4 

MBAND:  DATA 

INTEGER  ID  (NNOD)  . LM  (NDOF) 


C 

C--2.0  CALL  BCOND  TO  DO  THE  WORK 
C 

CALL  BCOND (IA (MPID) . NNOD, ERR) 
RETURN 


MBANEL  * 0 

DO  10  I -1. NDOF 

II  - ABS  (ID  (LM  (I)  ) ) 

DO  10  J-l.NDOF 
JJ  - ABS  (ID  (LM  ( J)  ) ) 

10  MBANEL  - MAX  (MBANEL.  ABS  (II-JJ+1) ) 

IF  (MBANEL  . GT . MBAN)  THEN 

WRITE  (NTM,  2000)  MBANEL,  MBAN 

WRITE  (NOT.  2000)  MBANEL.  MBAN 


ERR-. TRUE. 

ENDIF 


MAXBAN  - MAX (MAXBAN, MBANEL) 

RETURN 

2000  FORMAT  ( ' ***•  ERROR:  Bandwidth  exceeded  for  current  element.*/ 
. ' Bandwidth  required  :*,I5/ 

Bandwidth  available: ', 15) 


END 


C ADDCM 

SUBROUTINE  ADDCM  (ID, LM. ELM. SYS . NDOF. NSI ZE. NNOD. NEQN . MBAN) 

C-SUB:  ADDCM  - ADDS  ELEMENT  MATRIX  TO 
C COMPACT  FORM  OF  SYSTEM  MATRIX 


REAL*  8 ELM  (NSI ZE, NSIZE) , SYS  (NEQN, MBAN) 
INTEGER  ID  (NNOD)  .LM  (NSIZE) 

DO  20  I -1.  NDOF 
II  « ABS  (ID  (LM  (I))) 

DO  10  J-l.NDOF 

JJ  « ABS  (ID  (LM  ( J)  ) ) - II  +1 

IF ( JJ. LE. 0)  GO  TO  10 

SYS  (II , JJ)  - SYS  (II . JJ)  ♦ ELM(I.J) 

10  CONTINUE 
20  CONTINUE 

RETURN 

END 


C' 


2100  FORMAT (/'  — NOTE:  Boundary  condition  data  not  found.'/ 

All  nodes  assumed  to  be  flux-prescribed  nodes.'/ 
Equation  numbers  = node  numbers.') 


END 

C 

BCOND 

SUBROUTINE  BCOND (ID, NNOD. ERR) 

C-SUB: BCOND  - READS  AND  GENERATES  BOUNDARY  CONDITION  FLAGS 
C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR, NIN. NOT, ND1 , ND2 , ND3 , ND4 
C 

C BCOND:  DATA  4 COMMON  STORAGE 

C 

CHARACTER  ENDFLAG*3.  BC*1 
DIMENSION  ID  (NNOD).  IJK (3) 

LOGICAL* 1 ERR 

SAVE  /IOLIST/ 

C 

C— 1.0  WRITE  BOUNDARY  CONDITIONS  HEADER 
C 

WRITE (NOT.  2100) 

2100  FORMAT (/ 

. SX' Negative  Eqtn-#  - part  of  tenperat ure-prescrlbed  boundary.'/ 
,6x 'Positive  Eqtn-#  - part  of  flux-prescribed  boundary.*// 

.SX. 'Node  Eqtn-#  Node  Eqtn-#  Node  Eqtn-#  Wode  Eqtn-# 
Node  Eqtn-#  ' ) 

C 

C — 2.0  INTERPRET  LINE  OF  DATA 
C 

20  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH ('  ' .ENDFLAG. 3.1) 

IF ( ENDFLAG . EQ . ' END ' ) GO  TO  30 

C DO  IT 

CALL  FREEH  CC'.BC,  1,1) 

IF(BC.EQ. ’Q’ .OR.BC.EQ.'T')  THEN 
I SIGN  - 1 
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IF (BC.EQ. *T* ) ISIGN  « -1 
CALL  FREEI (’  MJK(1),3) 

IF (I JK(2) .EQ. 0)  IJK(2)-IJK(1) 

IF(IJK(3)  .EQ.O)  IJK(3)-1 
DO  24  N-I  JK  (1 ) , I JK  (2 ) , I JK  (3  ) 

IF(N.LE.O.OR.N.GT.NNOO)  THEN 
WRITE (NTH. 2200)  N 
WRITE (NOT, 2200)  N 
ERR =. TRUE. 

GO  TO  20 
ENDIF 

NTDOF  - ABS  (ID(N)  ) 

C SET  TDOF  TO  NEG.  VALUE  FOR  ALL  NOOES  CONSTRAINED  EQUAL 

C TO  AVOID  PROBLEMS  OF  MULTIPLE  BC  SPECIFICATION  USE  ABS () 

C SO  LAST  REDUNDANT  BC  SPECIFICATION  GOVERNS 


WRITE  (NTM, 2300) 

WRITE (NOT, 2300) 

RETURN 

ENDIF 

CALL  FREETY 
C 

C— 4.0  CALL  CONV  TO  DO  THE  WORK 
C 

CALL 

CONV  (IA(MPID)  , IA  (MPXYZ)  , IA  (MPK)  , I A (MPCH)  , NNCO,  NEQN  , MBAN  , ERR) 
RETURN 

2300  FORMAT (/’  — NOTE:  Convection  boundary  data  not  found.') 

END 


DO  24  NC-l.NNOD 
NCDOF-ABS (ID (NC) ) 

IF (NCDOF.EQ. NTDOF)  ID  (NC) «ISIGN*NCDOF 
24  CONTINUE 
ELSE 

WRITE  (NTM,  2210)  BC 
WRITE  (NOT,  2210)  BC 
ERR- . TRUE . 

ENDIF 
GO  TO  20 
C 

C— 3.0  REPORT  BOUNDARY  CONDITIONS  IF  NO  ERROR  ENCOUNTERED 
C 

30  IF (ERR)  RETURN 

WRITE  (NOT,  2300)  (N . ID (N) , N-l . NNOD) 


C CONV 

SUBROUTINE  CONV (ID. XYZ. K, CH. NNOD, NEQN . MBAN. ERR) 

C-SUB : CONV  - READS  AND  GENERATES  CONVECTION  BOUNDARY  DATA 
C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C 

C CAL - SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR. NIN, NOT. ND1 , ND2 , ND3 . ND4 
C 

C CONV:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 K (NEQN. MBAN) . CH (NNOD) . XYZ(NNOO,3) 

INTEGER  ID (NNOD) , LMNEW (2) . LMDLD (2) 

LOGICAL* 1 ERR 

CHARACTER  ENDFLAG*  3 . REG* 3 


RETURN 

2200  FORMAT (’  ••••  ERROR:  (Generated)  Node’, 15,’  Is  out  of  range.') 
2210  FORMAT  ('  *•••  ERROR:  Boundary  condition  ’.Al,'  not  available.') 
2300  FORMAT ((4X, 5(16, IX, 16, 2X) ) ) 

END 


C TREE:  CONV  - TREE  OF  SUBROUTINES  TO  ESTABLISH  CONVECTIVE  B.C.s 


SUBROUTINE  RCONV (MPID, MPXYZ . MPK , MPCH. NNOD. NEQN , MBAN . ERR ) 
C-SUB :RCONV  - ROOT  SUBROUTINE  TO  CONV 


C 

C 

C-; 

C 


ESTABLISHES  CONVECTIVE  BOUNDARY  CONDITIONS 


DICTIONARY  OF  VARIABLES 


VARIABLE  DESCRIPTION-- 

VARIABLES  PASSED  TO  SUBROUTINE 


C ERR 

C NNOD 

C NEQN 

C MBAN 

C MPCH 

C MPK 

C MPXYZ 

C MPID 

C CH (NNOD) 

EXCITATION 
C K (NEQN. MBAN) 

C ID (NNOD) 

C XYZ  (NNOD,  3) 

C 


DO-WHILE  TERMINATOR  FLAG 
NUMBER  OF  NODES  IN  SYSTEM 
NUMBER  OF  (SYSTEM)  EQUATIONS 
(HALF)  BANDWIDTH  OF  SYSTEM  EQUATIONS 
POINTER  TO  CH (NNOD)  IN  BLANK  COMMON 
POINTER  TO  K (NEQN)  IN  BLANK  COMMON 
POINTER  TO  XYZ  (NNOD, 3)  IN  BLANK  COMMON 
POINTER  TO  ID (NNOD)  IN  BLANK  COMMON 
COEFS.  FOR  COMPUTING  EFFECTIVE  CONVECTIVE 

SYSTEM  CONDUCTANCE  TRANSFER  MATRIX  (COMPACT  FORM) 
EQUATION  NUMBER/B.C.  FLAG  ARRAY; 

COORDINATE  ARRAY  (ORDERED  BY  NODE  NUMBER) 

ABS (ID  (N))  - EQUATION  NUMBER  OF  NODE  N 

SIGN (ID (N) ) - -1  * TEMPERATURE  PRESCRIBED  NODE 

SIGN (ID (N))  - ♦!  - FLUX  PRESCRIBED  NODE 


C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  MTOT.NP. IA (1) 

COMMON  /DBSYS/  NUMA. NEXT. IDIR. IP (3) 

COMMON  /IOLIST/  NTM, NTR, NIN, NOT, ND1 , ND2 . ND3 . ND4 


C ICOND:  DATA  4 COMMON  STORAGE 

C 

LOGICAL*!  ERR 


SAVE  /DBSYS/, /IOLIST/ 

C 

C— 0.0  WRITE  HEADER 
C 

WRITE (NTM, 2000) 

WRITE (NOT, 2000) 

2000  FORMAT  (/’  — - CONVECTIVE  BOUNDARY  CONDITIONS') 


C 

• C— 1.0  DEFINE  ARRAY 

C 

CALL  DELETE ('CH  ') 

CALL  DEFINE ('CH  MPCH. NNOD. 1 ) 
C 

c— 2.0  INITIALIZE  ARRAYS 
C 

CALL  ZEROR (IA (MPCH) .NNOD, 1) 

C 

C— 3.0  FIND  SEPARATOR  "CONVECT" 

C 

CALL  FINDNC  CONVECT ’,7.  KEY) 

IF (KEY . EQ. 1 ) THEN 


SAVE 

NDOF  = 2 


C 

C— 1.0  WRITE  CONVECTION  B.C.  HEADER 
C 

WRITE (NOT,  2100) 

2100  FORMAT (/ 

. €X ' D « Thickness  for  REG-PLN;  Subtended  angle  for  REG=AXI . * // 

.'  Surf  I J H-Coef  D REG') 

C 

C— 2.0  GET  FIRST  SURFACE  SEGMENT,  FORM  4 MODIFY  SYSTEM  K 
C 

CALL  FREE 

CALL  FREETY 

CALL  FREEI ('  ’ .NOLD, 1) 

CALL  FREEI ('I' ,LMOLD(l)  ,2) 

CALL  FREER ('H' .COEF, 1) 

D=1 . 0 

CALL  FREER ( 'D' .D.l) 

REG  * ' PLN ' 

CALL  FREEH  ('G'. REG. 3.1) 

CALL  CONVO ( ID, XYZ , NOLD, LMOLD . K. CH , COEF. REG. D. NNOD. NEQN . MBAN . ERR) 
C 

C--3.0  GET  NEXT  LINE  OF  DATA 
C 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH ( ' ' .ENDFLAG, 3.1) 

IF (ENDFLAG. EQ. 'END' ) THEN 
RETURN 
ENDIF 

C GET  NEW  ELEMENT  INFORMATION 

CALLFREEIC  '.NNEW.l) 

CALL  FREEI  ( ' I ' , LINEW  (1 ) ,2) 

INCR  * 1 

CALL  FREEI ('N' . INCR. 1) 

CALL  FREER ('H' , COEF, 1) 

0-1.0 

CALL  FREER ('D', D.l) 

REG  - 'PLN' 

CALL  FREEH  CG’.  REG.  3.1) 

C CHECK  NUMERICAL  ORDER 

IF (NNEW . L£ . NOLD ) THEN 
WRITE (NTM, 2300)  WNEM 
WRITE (NOT, 2300)  NNEW 
ERR-. TRUE. 

RETURN 

ENDIF 

C GENERATE  MISSING  SURFACE  SEGMENTS 

DO  34  N— NOLD+1 . NNEW- L,  1 
DO  32  I -1. NDOF 

32  LMOLD  (I)  - LMOLD  (I)  ♦INCR 

34  CALL  CONVO (ID, XYZ. N, LMOLD, K.CH, COEF. REG, D, NNOD. NEQN, MBAN. ERR) 

C DO  NEW  SURFACE  SEGMENT 

HOLD  - NNEW 
DO  36  I-l.NDOF 
36  LMOLD  (I)  - LMJEM  (I ) 

CALL  CONVO (ID. XYZ. NOLD, LMOLD. K. CH.COEF. REG, D. NNOD. NEQN. MBAN.  ERR) 
GO  TO  30 

2300  FORMAT (’  ••••  ERROR : Surface  segment  number ',15,'  out  of 
order. ') 


D -1  1 


RLSOLV 


END 


C CONVO 

SUBROUTINE  CONVO (ID. XYZ . NSEG. LM. K.CH.COEF. REG. D.  NNOD.  NEQN. 

■►MEAN,  ERR) 

C-SUB : CONVO  - REPORTS  CONVECTIVE  SURFACE  SEGMENT  INFORMATION, 

C FORMS  COEFS,  CH (NNOD) . FOR  COMPUTING  EFFECT . CONV . EXCIT . 

C MODIFIES  K MATRIX  FOR  CONVECTION  BOUNDARY  CONDITIONS 

C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM, NTR, NIN , NOT, ND1.ND2.ND3.ND4 
C 

C rORM:  DATA  4 COWON  STORAGE 

C 

REAL*  8 K (NEQN , MEAN)  , CH (NNOD) , XYZ (NNOD. 3) . S(2,2) 

INTEGER  ID  (NNOD) ( LM(2) 

LOGICAL*l  ERR 
CHARACTER  REG* 3 

SAVE 

C 

C— 1.0  REPORT  ELEMENT  INFORMATION  TO  OUTPUT  DATA  FILE 
C 

WRITE (NOT. 2000)  NSEG. LM (1) ,LM(2) .COEF.D.REG 
C 

C — 2.0  COMPUTE  COEFFICIENT. CH.  FOR  COMPUTING  EFFECTIVE  CONVECTIVE  FLUX 
C 

C ERROR  TRAP  FOR  REGION  TYPE 

IF (REG. EQ. * PLN ' .OR . REG. EQ. ' AXI ' ) THEN 
CONTINUE 
ELSE 

WRITE  (NTM.  2200)  REG 
WRITE  (NOT.  2200)  REG 
ERR= . TRUE . 

RETURN 
END  IF 

C ERROR  TRAP  FOR  OUT-OF-RANGE  NODES 

DO  200  N=1 , 2 
NN  * LM  (N) 

IF (NN . LE. 0 . OR . NN . GT .NNOD)  THEN 
WRITE  (NTM.  2210)  NN 
WRITE  (NOT,  2210)  NN 
ERR*. TRUE. 

RETURN 

C ERROR  TRAP  FOR  TEMP -PRESCRI BED  NODES 

ELSEIF (ID (NN) .LE. 0)  THEN 
WRITE  (NTM.  2220)  NN 

WRITE  (NOT.  2220)  NN 


ERR*. TRUE. 


RETURN 

ENDIF 

200  CONTINUE 
N1  * LM  (1 ) 

N2  * LM  (2) 

C COMPUTE  SURFFACE  SEGMENT  LENGTH 

XL  * SORT ( (XYZ (Nl,l) -XYZ (N2,l))**2  ♦ 

. (XYZ  (Nl. 2) -XYZ  (N2.2)  )**2) 

C COMPUTE  CH 

IF (REG.EQ. 'PLN' ) THEN 

CH (Nl ) * CH  (Nl)  ♦ D*COEF*XL/2. 0 
CH(N2)  « CH(N2)  ♦ D*COEF*XL/2 . 0 
ELSEIF (REG.EQ. ’AXI ’ ) THEN 

CH (Nl ) * CH(N1)4 (D*COEF*XL/6. 0)*  (2.0*XYZ  (Nl. 1)  + XYZ  (N2. 1 ) ) 
CH(N2)  * CH(N2)4*  (D*COEF*XL/6 . 0)*  (XYZ  (Nl.  1)  ♦ 2 . 0*XYZ  (N2. 1 ) ) 
ENDIF 
C 

C — 3.0  MODIFY  CONDUCTANCE  TRANSFER  MATRIX.  K 
C 

C COMPUTE  CONVECTIVE  CONTRIBUTION  TO  K 

IF (REG.EQ. ’PLN' ) THEN 
S(l,l)  - D*COEF*XL/3 . 0 
S (1 . 2)  - D*COEF*XL/4. 0 
S (2 , 1 ) - S (1 , 2) 

S (2 , 2)  - S(l.l) 

ELSEIF (REG.EQ. ’AXI’)  THEN 

S(l.l)  - (D*COEF*XL/4.0)* (XYZ(Nl.l)  ♦ XYZ  (N2 . 1 ) /3 . 0) 

S (1 , 2)  - (D*COEF*X1/12.0)*(XYZ(N1,1)  ♦ XYZ(N2,1)) 

S (2 , 1 ) - S(1.2) 

S (2 , 2)  - (D*COEF*XL/ 4.0)* (XYZ (Nl , 1 ) /3  ♦ XYZ(N2,1>) 

ENDIF 

C ADD  CONVECTIVE  CONTRIBUTION  TO  SYSTEM  ARRAYS 

CALL  ADDCM (ID, LM, S , K, 2, 2, NNOD, NEGN . MBAN) 

RETURN 

2000  FORMAT (6X,I4,2I5.5X. 2G10. 3, 3X, A3) 

2200  FORMAT ( ’ ***•  ERROR:  Region  ' , A3, * is  not  defined.’) 

2210  FORMAT ('  ••••  ERROR:  (Generated)  Node’,15, ’ la  out  of  range.’) 

2220  FORMAT (’  •*••  ERROR:  (Generated)  Node’. 15, 

.*  la  not  part  of  flux-prescribed  boundary') 


SUBROUTINE  RLSOLV (MPID. MPK, MPC, MPEG. MPCH . NNOD. NEQN . MBAN  . ERR) 
C--SUB : RLSOLV  - ROOT  SUBROUTINE  TO  ALL  LINEAR  SOLVER  SUBROUTINES 
C 

C DICTIONARY  OF  VARIABLES 

C 

C VARIABLE 

C 

C ERR 

C NNOD 

C NEQN 

C MBAN 

C IWRT 

C IPRT 

C NPRT 

C MPNPR 

C MPCH 

C MPEO 

C MPC 

C MPK 

C MPID 

C NPR(NPRT) 

C CH  (NNOD) 

EXCITATION 
C EO(NNOO) 

C C (NEQN . MBAN) 

C K (NEQN , MBAN) 

C ID (NNOD) 

C 

c 
c 

c 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  MTOT.NP, IA(1) 

COMMON  /DBSYS/  NUMA, NEXT. IDIR, IP  (3) 

COMMON  /IOLIST/  NTM. NTR. NIN , NOT. ND1 . ND2 , ND3 . ND4 
C 

C RFORM:  DATA  4 COWON  STORAGE 

C 

LOGICAL* 1 ERR,  NOSOLV 
INTEGER  ITIME1. ITIME2 

SAVE  /DBSYS/, /IOLIST/ 

C 

C--0.0  WRITE  SOLUTION  HEADER  4 WRITE  SYSTEM  MATRICES  TO  <f llename> . SYS 
C 

WRITE  (NTM. 2100) 

WRITE (NOT. 2100) 

2100  FORMAT  (/,’  -=  SOLUTION’) 

CALL  SYSAVE (IA(MPK) , IA (MPC) .NEQN. MBAN) 

NOSOLV  * .TRUE. 

C 

C — 1.0  DEFINE  E(NECN) 

C 

CALL  DELETECE  ’) 

CALL  DEFINECE  ’ . MPE . NEQN.  1 ) 

C 

C— 2.0  GET  REPORT  CONTROL  INFORMATION 
C 

CALL  REPORT (MPNPR, NPRT. IPRT. IWRT, NNOD. ERR) 

IF (ERR)  RETURN 
C 

C--3.0  LOOK  FOR  SEPARATOR  "STEADY" 

C 

CALL  FINDN (’STEADY’ , 6. KEY) 

IF (KEY . EQ. 0)  THEN 
NOSOLV*. FALSE. 

CALL  ZEROR (IA (MPE) .NEQN, 1) 

CALL  TIME (ITIME1 ) 

CALL  STEADY (IA (MPID) . IA (MPK) . IA (MPE) , IA (MPEO) .IA(MPCH) , 

♦ IA (MPNPR) .NPRT. IWRT. NNOD, NEQN. MBAN. ERR) 

CALL  TIME (ITIME2) 

WRITE (NTM. 2300)  (ITIME2-ITIME1 ) 

WRITE (NOT. 2300)  (ITIME2-ITIME1 ) 

ENDIF 

2300  FORMAT (/’  — NOTE:  Solution  time:  ',15, • aeconda.’) 


DESCRIPTION 

DO-WHILE  TERMINATOR  FLAG 
NUMBER  OF  NODES  IN  SYSTEM 
NUMBER  OF  (SYSTEM)  EQUATIONS 
(HALF)  BANDWIDTH  OF  SYSTEM  EQUATIONS 
WRI TE- SOLUTION -TO- FILE  FLAG  (l*DO  IT) 

OUTPUT  PRINT  INTERVAL 

NUMBER  OF  NOOES  SELECTED  FOR  SOLUTION  OUTPUT  PRINT 
POINTER  TO  NPR  (NPRT)  IN  BLANK  COMMON 
POINTER  TO  CH(NNOO)  IN  BLANK  COWON 
POINTER  TO  RO  (NNOD)  IN  BLANK  COWON 
POINTER  TO  C (NECN . MBAN)  IN  BLANK  COWON 
POINTER  TO  K (NEQN, MBAN)  IN  BLANK  COMMON 
POINTER  TO  ID  (NNOD)  IN  BLANK  COWON 
LIST  OF  NODES  SELECTED  FOR  SOLUTION  OUTPUT  PRINT 
COEFS.  FOR  COMPUTING  EFFECTIVE  CONVECTIVE 

INITIAL  EXCITATION  VECTOR  (STORED  BY  NODE) 

SYSTEM  CAPACITY  MATRIX  (COMPACT  FORM) 

SYSTEM  CONDUCTANCE  TRANSFER  MATRIX  (COMPACT  FORM) 
EQUATION  NUMBER /B.C.  FLAG  ARRAY; 

ABS  (ID  (N) ) « EQUATION  NUMBER  OF  NODE  N 

SIGN (ID (N))  - -1  - TEMPERATURE  PRESCRIBED  NOOE 

SIGN (ID  (N) ) - ♦!  - FLUX  PRESCRIBED  NODE 


C 

C— 4.0  LOOK  FOR  SEPARATOR  "HARMON" 
C 

CALL  riNDNC  HARMON’,  6,  KEY) 

IF  ( KEY . EQ . 0 ) THEN 
NOSOLV*. FALSE. 

CALL  ZEROR (IA (MPE) .NEON, 1) 

CALL  TIME (ITIME1 ) 


END 


C»* *•*••*••»*•«•••**»•*•••••••**•••*••••*•*•**•***•• •*••***»••*•< 

C TREE:  LSOLV  - TREE  OF  SUBROUTINES  TO  SOLVE  LINEAR  SYSTEM  EQTNS 

c 


CALL  HARMON  (IA  (MPID)  . IA  (MPK)  . IA  (MPC)  . IA  (MPE)  . IA  (MPEO)  . IA  (MPCH)  . 
. IA (MPNPR) . IA (MPK) .IA(MPE) . NPRT. IWRT. NNOD, NEQN . MBAN . ERR) 
CALL  TIME (ITIME2 ) 

WRITE  (NTM. 2300)  (ITIME2-ITIME1 ) 


-1  2 


D 


WRITE (NOT. 2300)  (ITIME2-ITIME1) 

END  IF 
C 

c — 5.0  LOOK  FOR  SEPARATOR  "PREDICT" 

C 

CALL  FINDN (’PREDICT*.  7. KEY) 

IF  (KEY . EQ. 0)  THEN 
NOSOL V*. FALSE. 

CALL  ZEROR(IA(MPE) .NECN.l) 

CALL  TIKE (ITINE1 ) 

CALL  PREDIC (IA (MPID) . IA  (MPK) . IA  (HPC) . IA (MPE) . IA (KPEO) . IA (MPCH) , 
. IA(KPNPR) .NPRT.IPRT, IWRT. NNOD, NEON . MBAN. ERR) 

CALL  TIKE (ITIME2 ) 

WRITE  (NTH. 2300)  (ITIME2-ITIKE1) 

WRITE (NOT. 2300)  (ITIHE2-ITIKE1) 

END  IF 
C 

C— 6.0  REPORT  NO  SOLUTION  REQUEST  ... 

C 

IF (NOSOLV)  THEN 
WRITE  (NTH.  2600) 

WRITE  (NOT.  2600) 

ENDIF 

CALL  FCLOSE  (ND1 ) 


RETURN 

2600  FORMAT ( / ' — NOTE:  No  solution  control  data  found.') 

END 

C SYSAVE 


SUBROUTINE  SYSAVE (K. C. NEQN . MBAN ) 

C-SUB : SYSAVE  - WRITES  (K)  4 (C)  TO  <f llenarae> . SYS 
C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLI ST/NTM, NTR . NIN . NOT . ND1 , ND2 , ND3 . ND4 
C 

C SYSAVE:  DATA 

C 

REAL*  6 K (NEQN , MEAN ) . C (NEQN . MBAN ) 

C 

CALL  FOPEN (ND1. ’SYS’, ’NEW’ ) 

WRITE (ND1)  ((K(N.M) , N=1 , NEQN) . M=1 , MBAN ) , 

+ ((C(N,M),N*1,  NEQN ) , M=1 . MBAN) 

RETURN 

END 


c STEADY 

SUBROUTINE  STEADY (ID. K, E, EO. CH, NPR. NPRT, IWRT. NNOD. NEQN , MBAN . ERR) 
C-SUB:  STEADY  - FORMS  (E)  OF  (K)(T)  * (E) 

C CALL  SOLVER  AND  REPORTS  RESULTS 

C SOLUTION  (T)  WRITTEN  OVER  (E) 


IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL -SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLI ST/NTM. NTR. NIN , NOT, ND1 , ND2 , ND3 , ND4 
C 

C STEADY:  DATA  4 COMMON  STORAGE 

C 

REAL* 8 K (NEQN. MBAN) . E (NEQN) , EO (NNOD) , CH (NNOD) 
INTEGER  ID (NNOD).  NPR  (NPRT) , IJK(3) 

LOGICAL *1  ERR.  TDOF 
CHARACTER  ENDFLAG*  3 

SAVE  /IOLIST/ 

C 

C— 1.0  WRITE  HEADER 
C 

WRITE (NTH, 2100) 

WRITE (NOT, 2100) 

2100  FORMAT (/’  STEADY  STATE  SOLUTION’) 

CALL  FREETY 
C 

C— 2.0  SCALE  (K)  FOR  TEMP -PRESCRIBED  NODES 
C 

DO  20  I -1. NEQN 

20  IF (TDOF (I. ID. NNOD))  K(I,1)  - K (I . 1) *1 . OE+15 
C 

C~ 3.0  GET  T,  Q.  AND/OR  TEX  4 FORM  { E ) 

C 

WRITE (NTM. 2300) 

WRITE (NOT. 2300) 

WRITE (NOT,  2302) 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END" 

CALL  FREEH (’  ’. ENDFLAG. 3. 1 ) 

IF (ENDFLAG. EQ. ’END’ ) GO  TO  40 

C INTERPRET  LINE  AND  GENERATE  DATA 

Q - 0.0 

CALL  FREER ('O’ .Q.  1) 

T = 0.0 

CALL  FREER (’T* . T, 1) 

TEX  - 0.0 


CALL  FREER  CX*  .TEX.  1) 

CALL  FREEI ( ' ’.IJK(1),3) 

IF  (I  JK  (2)  .EQ.  0)  I JK  (2)  “I  JK  (1 ) 

IF  (I  JK  (3)  . EQ.  0)  I JK  (3)  =1 
DO  32  N«IJK(1)  ,IJK(2).IJK(3) 

IF (N . LE . 0 . OR . N . GT . NNOD ) THEN 
WRITE  (NTM. 2310)  N 
WRITE (NOT. 2310)  N 
ERR* . TRUE . 

GO  TO  30 
ENDIF 
NN  - ID  (N) 

NNN  - ABS  (NN) 

C TEMPERATURE-PRESCRIBED  NODES 

IF(NN.LT.O)  THEN 

WRITE (NOT. 2320)  N.T 
E (NNN)  - T*K (NNN , 1 ) 

C FLUX -PRESCRIBED  NODES 

ELSEIF  (NN.GT.0)  THEN 

WRITE (NOT.  2330)  N. EO (N) . Q. TEX 
E (NNN)  ■ E (NNN)  ♦ EO(N)  + Q + TEX*CH(N) 

ENDIF 

32  CONTINUE 
GO  TO  30 
C 

C — 4 . 0 SOLVE 
C 

40  IF (ERR)  RETURN 

CALL  SYMBC(K.E.NECN.  MBAN.  NEQN.  1) 

C 

C— 5.0  REPORT  SOLUTION  (T) 

C 

C WRITE  TO  <FILENAME> . TMP  OPTION 

IF (IWRT . EQ. 1)  THEN 
TIME  * 0.0 

CALL  FOPEN (ND1. ’TMP  ’.’NEW’) 

WRITE (ND1)  TIME,  (E  (N)  . N-l . NEQN) 

CALL  FCLOSE (ND1) 

ENDIF 

C PRINTABLE  OUTPUT 

WRITE (NTM. 2500) 

WRITE  (NTM. 2510)  (NPR (N)  . E (ABS (ID (NPR (N) ) ) ) . N»1.NPRT) 

WRITE (NOT. 2500) 

WRITE  (NOT.  2510)  (NPR  (N)  . E (ABS  (ID  (NPR  (N)  ) ) ) ,N=1.NPRT) 

RETURN 

2300  FORMAT (/’  — EXCITATION’) 

2302  FORMAT ( / 

.6X'Node  Terap.  Int.Flux  Ext. Flux 

. Ext . Temp. ' ) 

2310  FORMAT ( ’ •*•*  ERROR:  (Generated)  Node'. 15,'  Is  out  of  range.') 
2320  FORMAT ( 6X , I 4 , 9X , G1 0 . 3 ) 

2330  FORMAT (6X, 14, 19X, 3 (5X.G10.3) ) 

2500  FORMAT (/ ’ — RESPONSE'// 

.6X,'Node  Tenp . Node  Terap.  Node  Terap.  Node 

Terap.  Node  Tern?.') 

2510  FORMAT ( (6X. 5 (I4.G11 .3) ) ) 

END 


C HARMON 

SUBROUTINE  HARMON (ID. K.C.E. EO. CH. NPR, KSTAR. ESTAR. 

.NPRT. IWRT, NNOD. NEQN. MBAN. ERR) 

C-SUB:  HARMON  - FORMS  (K*)(T*>  « ( E* ) 

C CALL  SOLVER  AND  REPORTS  RESULTS 

C COMPLEX  WRITTEN  OVER  DOUBLE  PRECISION  STORAGE  CELLS: 

C (K* 1 WRITTEN  OVER  [K] 

C (E*)  WRITTEN  OVER  (E) 

C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/NTM.NTR.NIN.NOT.NDl.ND2.’ND3.ND4 

c 

C HARMON:  DATA  4 COMMON  STORAGE 

C 


REAL* 8 K (NEQN. MBAN) , C (NEQN, MBAN) , E (NEQN) , EO (NNOD) ,CH (NWOO)  . 
+0(2) ,T (2) , TEX (2) 

INTEGER  ID  (NNOD)  , NPR  (NPRT)  . IJK(3) 

COMPLEX  KSTAR  (NEQN. MBAN ) , ESTAR (NEQN) . QSTAR.  TSTAR,  TXSTAR , 1 
LOG  I CAL*  1 IRR,  TDOF 
CHARACTER  ENDFLAG* 3 
REAL  ZLAG 

SAVE  /IOLIST/ 

DATA  PI/3.141592654/ 

C 

C— 1.0  WRITE  HEADER 

C 

WRITE (NTM, 2100) 

WRITE (NOT.  2100) 

2100  FORMAT (/’  STEADY  HARMONIC  SOLUTION’) 

CALL  FREETY 
C 

C--2.0  GET  FREQUENCY 
C 


D -1  3 


CALL  FREE 
CALL  FREETY 
W - 0.0 


2610  FORMAT ( (6X,3  (14, 2G 10.3)) ) 


CALL  FREER ('F' ,W. 1) 

WRITE (NOT, 2200)  W 

2200  FORMAT (/ ' Frequency  of  excltat ion: * ,G10 .3, ' cycles /time ' ) 

M - W*  2 . 0*PI 
C 

C— 3.0  GET  (K)  4 (C),  FORM  (K*),  6 SCALE  (K* ) FOR  TEMP -PRESCRIBED 

NODES 

C 

REWIND  ND1 

READ (ND1 ) < (K(N,M) .N-l. NEQN)  ,M-1, MBAN)  , 

. (<C(N,M)  .N-l.NECN)  , M=1 , MBAN ) 

DO  30  I-l.NEQN 
DO  30  J-l.MBAN 

30  KSTAR  (I , J)  - CMP  LX  (K  (I , J)  , W'C(I,J|) 

DO  32  I-l.NEQN 

32  IF (TDOF  (I, ID.NNOD) ) KSTAR(I.l)  - KSTAR (1,1)* (1 . 0E+ 1 5 , 0.0) 

C 

C — 4.0  GET  T,  Q.  AND/OR  TEX  4 FORM  (E) 

C 

WRITE (NTM, 2400) 

WRITE (NOT, 2400) 

WRITE (NOT, 2402) 

40  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END” 

CALL  FREEH ('  ', ENDFLAG . 3.  1 ) 

IF (ENDFLAG . EQ. ' END ' ) GO  TO  50 
C INTERPRET  LINE  AND  GENERATE  DATA 

C CONVERT:  Z - MAG* COS (-LAG*W)  * 1 MAG*SIN (-LAG*W) 

Q (1 ) - 0.0 
Q(2)  - 0.0 
CALL  FREER  <’Q’ ,Q.  2) 

QSTAR  - CMPLX (Q(l) *COS  (-Q(2) *W) ,Q(1) *SIN  (-0(2) *W) ) 

T (1 ) - 0.0 
T (2 ) - 0.0 
CALL  FREER ( 'T* , T,  2) 

TSTAR  = CMPLX (T (1) *COS (-T (2) *W) ,T(1) *SIN (-T (2) *W) ) 

TEX(l)  =0.0 
TEX (2)  - 0.0 
CALL  FREER (’X’, TEX. 2) 

TXSTAR  = CMPLX (TEX (1) *COS (-TEX(2)*W), TEX (1 ) *SIN (-TEX (2) *W) ) 

CALL  FREEI  ( * ',IJK(1),3) 

IF (IJK (2) .EQ. 0)  IJK (2) =IJK (1) 

IF  (I JK  (3)  . EQ.  0)  IJK  (3)  =1 
DO  42  N=IJK (1) , IJK (2) , IJK (3) 

IF (N.LE.O.OR.N.GT.NNOD)  THEN 
WRITE  (NTM.  2410)  N 
WRITE  (NOT,  2410)  N 
ERR= . TRUE . 

GO  TO  40 
ENDIF 
NN  - ID  (N) 

NNN  = ABS (NN) 

C TEMPERATURE-PRESCRIBED  NODES 

IF(NN.LT.O)  THEN 

WRITE  (NOT.  2420)  N,T(1),T(2) 

ESTAR  (NNN ) - TSTAR* KSTAR (NNN . 1 ) 

C FLUX-PRESCRIBED  NODES 

ELSEIF (NN.GT. 0)  THEN 

WRITE  (NOT.  2430)  N.  EO  (N)  , Q (1)  . Q (2)  .TEX(l)  ,TEX(2) 

ESTAR  (NNN)  - ESTAR (NNN)  ♦ CMPLX (EO(N))  ♦ QSTAR  ♦ CH(N)*TXSTAR 
ENDIF 

42  CONTINUE 
GO  TO  40 
C 

C— 5.0  SOLVE 
C 

50  IF (ERR)  RETURN 

CALL  SYMBCX (KSTAR. ESTAR, NEQN. MBAN. NEQN. 1 ) 

C 

C — 6.0  REPORT  SOLUTION  { T ) 

C 

C WRITE  TO  <FILENAME> . TMP  OPTION 

IF  (IWRT.EQ.l)  THEN 
TIME  - 0.0 

CALL  FOPEN(NDl, 'TMP  ’,,NEW’) 

WRITE  (ND1)  TIME.  (ESTAR (N) , N=1 . NEQN) 

CALL  FCLOSE(NDl) 

ENDIF 

C PRINTABLE  OUTPUT 

WRITE (NOT.  2600) 

WRITE  (NOT.  2610)  (NPR  (N)  , ABS  (ESTAR  (ABS  (ID  (NPR  (N)  ) ) ) ) , 

. SLAG (ESTAR (ABS (ID (NPR (N) ) )) , W) . N-l.NPRT) 

WRITE  (NTM,  2600) 

WRITE  (NTM,  2610)  (NPR  (N)  , ABS  (ESTAR  (ABS  (ID  (NPR  (N) ))))  , 

. ZLAG (ESTAR (ABS (ID  (NPR  (N) ) ) ) . W) , N=1 . NPRT) 

RETURN 

2400  FORMAT (/’  — EXCITATION’) 

2402  FORMAT  (/ 

.23X, f Tenperature  Ext.  Flux  Ext.  Tenp*/ 

. 6X, 'Node  Int.Flux  Arap.  Lag  Artp.  Lag 

Arap . Lag ' ) 

2410  FORMAT ( ' **•*  ERROR:  (Generated)  Node', 15. * la  out  of  range.') 
2420  FORMAT (6X, 14, 10X, 2G10 . 3) 

2430  FORMAT ( 6X . I 4 . G1 0 . 3 , 20X , 4G1 0.3) 

2600  FORMAT (/ ' — RESPONSE’// 

.13X’ Tenperature  Temperature 

Tenperature’ / 

. 6X’Node  Anp.  Lag  Node  Anp.  Lag 

. Node  Anp.  Lag') 


END 

C ZLAG 

FUNCTION  ZLAG(Z.W) 

C-FUN : ZLAG  - COMPUTES  LAG-ARG/FREQ  OF  COMPLEX  NUMBER 
COMPLEX  Z 
REAL* 8 W 

REAL  ZLAG, PI, IM, RE 
DATA  PI/3.141592654/ 

IM  - AIMAG(Z) 

RE  - REAL (Z) 


IF( (ABS (RE) .LT.1.0E-6) .OR. (W . LT . 1 . 0D-6) ) THEN 
ZLAG  - 0.0 
ELSE 

ZLAG  - -ATAN (IM/RE) 

IF (IM . GT. 0.0. AND . RE .GT .0.0)  ZLAG  - 2.0*PI+ZLAG 
IF (RE . LT .0.0)  ZLAG  - PI+ZLAG 
ZLAG  - ZLAG /REAL (W) 

ENDIF 

RETURN 

END 


C 

PREDIC 


SUBROUTINE  PREDIC (ID, K. C. E, EO.CH. NPR. NPRT. IPRT, IWRT. 
♦NNOD, NEQN,  MBAN,  ERR) 

C-SUB:  PREDIC  - PREDICTOR -CORRECTOR  1ST  O.D.E.  EQUATION  SOLVER 
C BASED  ON  SOLVER  IN  *HEAT*  BY  R.L. TAYLOR  - U.C. 

BERKELEY 

C SOLVES  EQUATION; 


(K)(T)  ♦ (C) (dT/dt ) - ( E (t ) ) 

FOR  GENERAL  EXCITATION,  { E (t ) > . DEFINED  BY  A SET  OF  PIECEWISE 
LINEAR  FUNCTIONS. 

METHOD  BASED  ON  DIFFERENCE  APPROXIMATION; 

(T) r+l  - (T) n ♦ (1-a) DT{dT/dt ) n ♦ (a) DT(dT/dt }n+l 


"alpha",  an  Integration  parameter 

0 corresponds  to  Forward  Difference  method 

1 corresponds  to  Backward  Difference  method 
1/2  corresponds  to  Crank-Nichol son  method 


(unstable) 

C DT  = time  step  increment 

C 

C DICTIONARY  OF  VARIABLES 

C INTERNAL  TO  SUBROUTINE 

C 

C VARIABLE  DESCRIPTION 

C NEFN 

NUMBER  OF  EXCITATION  FUNCTIONS  DEFINED 
C MPIEFN 

IEFN (NNOD, 2) 

: EXCITATION  FUNCTION  NUMBER  ARRAY 
C MPEFN 

EFN (NEFN.  2) 

: EXCITATION  FUNCTION  DATA 
C E (NEQN) 


: CURRENT  (E)  (ORDERED  BY  EQTN  #) 
C MPT 

T (NEQN) 


: CURRENT  (T)  (ORDERED  BY  EQTN  #) 

C MPTD 

TD  (NEQN) 

: CURRENT  (dT/dt)  (ORDERED  BY  EQTN  #) 

C MPTDD 

TDD (NEQN) 

: INITIAL  (d/dt (dT/dt ) ) TO  EST  TIME  STEP 

C IEFN (NNOD, 2)  LIST  OF  FUNCTION  NUMBERS  (ORDERED  BY  NODE  NUMBER) 

C IEFN (I. 1)  = TEMP  OR  EXT. FLUX  FUNCTION  NUMBER 

C IEFN (I, 2)  - CONVECTIVE  FLUID  TEMP.  FUNCTION  NUMBER 

C EFN (NEFN. 2)  EXCITATION  FUNCTION  DATA 

C EFN (1,1)  « OLD  VALUE 

C EFN (I, 2)  - NEW  VALUE 

C TOLD.  TNEW  EXCITATION  FUNCTION  DATA  TIMES 

C TIME  (1)  START  TIME 

C TIME (2)  END  TIME 

C TIME (3)  TIME  INCREMENT 


IMPLICIT  REAL*8  (A-H.O-Z) 

C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  MTOT.NP, IA(1) 

COMMON  /DBSYS/  NUKA. NEXT. IDIR, IP (3) 

COMMON  /IOLIST/  NTM. NTR . NIN . NOT . ND1 . ND2 . ND3 . ND4 

C 

C PREDIC:  DATA  4 COWON  STORAGE 

C 

REAL* 8 K (NEQN, MBAN) . C (NEQN, MBAN) . E (NEQN) , EO (NNOD) , CH (NNOD) 
INTEGER  ID (NNOD),  NPR (NPRT) 

CHARACTER  ENDFLAG* 3 
LOGICAL* 1 ERR 

COMMON  /LIST/  TIME (3) .ALPHA. TOLD. TNEW 
SAVE  /DBSYS/, /IOLIST/. /LIST/ 
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D 


C WRITE  HEADER 


CALL  FREETY 


C 

WRITE (NOT.  2000) 

2000  FORMAT (/*  — DYNAMIC  SOLUTION:  PREDICTOR-CORRECTOR 

INTEGRATION’ ) 

C 

C — 1.0  GET  SOLUTION  CONTROL  INFORMATION 
C 

CALL  FREE 
CALL  FREETY 

CALL  FREER ('E'. TIME (1 ) .3) 

ALPHA  « 0.75 

CALL  FREER (’A* .ALPHA.  1) 

IF(ALPHA. LT. 0.0. OR. ALPHA. GT. 1.0)  THEN 
WRITE  (NTH.  2100) 

WRITE  (NOT.  2100) 

ERR-. TRUE. 

RETURN 

2100  FORMAT (’  ERROR:  Alpha  «jst  be  In  range  0.0  to  1.0.') 

ENDIF 
NEFN  - 0 

CALL  FREEI (’N* .NEFN.l) 

IF (NEFN . LE . 0)  THEN 
WRITE  (NTM.2110) 

WRITE  (NOT.  2110) 

ERR- . TRUE . 

RETURN 

2110  FORMAT ( ' ****  ERROR:  One  or  more  excitation  functions 
. must  be  defined.') 

ENDIF 

C REPORT  CONTROL  INFORMATION 

WRITE (NOT, 2120)  TIME. ALPHA . NEFN 
2120  FORMAT (/’  — SOLUTION  CONTROL  INFORMATION’/ 


.’  Start  time  ’.G10.3/ 

End  time  ' .G10.3/ 

.'  Time  step  Increment  '.G10.3/ 

. ’ Integration  parameter,  alpha  '.G10.3/ 

. ' Number  of  excitation  functions  ...  ’,16) 


C 

C— 2.0  DEFINE  ARRAYS 
C 

CALL  DELETE ('IEFN') 

CALL  DELETE ( ’ EFN  ') 

CALL  DELETE (’T  ’) 

CALL  DELETE ( ' TD  ’) 

CALL  DELETE (’TDD  ’) 

CALL  DEFINE ('TDD  ' , MPTDD, NEQN , 1 ) 

CALL  DEFINE ( ’ TD  * , MPTD, NEQN, 1 ) 

CALL  DEFINE ('T  ’. MPT. NEQN . 1 ) 

CALL  DEFINE (’EFN  ’ . MPEFN, NEFN , 2 ) 

CALL  DEFINI (’IEFN’ . MPIEFN. NNOD, 2) 

NS TOR  « (IDIR-NEXT-20) 

IF (NSTOR.LT. 0)  THEN 

WRITE (NTM. 2200)  -NSTOR 

WRITE  (NOT.  2200)  -NSTOR 

ERR= . TRUE . 

RETURN 

2200  FORMAT ('  •***  ERROR:  Insufficient  storage  to  form  arrays  for 
. solution.'/ 

. ' Increase  IA (MTOT)  by  ’,16) 

ENDIF 

C INITIALIZE 

CALL  ZEROI (IA (MPIEFN) .NNOD. 2) 

CALL  ZEROR (IA (MPEFN) .NEFN. 2) 

CALL  ZEROR  (IA (MPT) .NEQN. 1) 

CALL  ZEROR (IA (MPTD) .NEQN. 1) 

CALL  ZEROR (IA (MPTDD)  , NEQN , 1 ) 

C 

C— 3.0  GET  EXCITATION  CONTROL  INFORMATION 
C 

CALL  ECNTRL  (ID.  IA  (MPIEFN)  . NNGD , ERR) 

IF (ERR)  RETURN 
C 

C— 4.0  GET  INITIAL  CONDITIONS 
C 

C 4.1  WRITE  HEADER 

C 

WRITE (NTM, 2400) 

WRITE (NOT,  2400) 

2400  FORMAT  (/’  — INITIAL  CONDITIONS') 

C 

C 4.2  FIND  SEPARATOR  "INITIAL" 

C 

CALL  FINDN  (’INITIAL '.7,  KEY) 

IF (KEY . EQ. 1)  THEN 
WRITE  (NTM.  2410) 

WRITE  (NOT.  2410) 

2410  FORMAT (/’  — NOTE:  Initial  condition  data  not  found.'/ 

.’  Initial  temperatures  assumed  to  be  zero.’) 

GO  TO  50 
ENDIF 


C 4.3  CALL  ICOND  TO  DO  THE  WORK 

C 

CALL  ICOND (ID, IA (MPT) .NNOD. NEQN, ERR) 

C 

C— 5.0  4 READ  (K)  4 (C)  FROM  <f  llenarae> . SYS 
C 

50  REWIND  ND1 

READ  (ND1 ) ( (K  (N,  M)  . N— 1,  NEQN)  . M-l , MB  AN)  , 

. ((C(N.M)  ,N-1  .NEQN)  .M-l.MBAN) 

C 

C--6.0  INITIALIZE  EXCITATION  FUNCTION  VALUES 
C 

CALL  FINDN ('EXCIT' . 5, KEY) 

IF (KEY . EQ. 1)  THEN 
WRITE (NTM. 2400) 

WRITE (NOT. 2600) 

2600  FORMAT ('  •••*  ERROR:  Separator  "EXCIT"  not  found.’) 

ERR- . TRUE . 

RETURN 

ENDIF 

DO  60  1-1.2 

CALL  GETEFN (IA (MPEFN) .TOLD. TNEW. NEFN. ENDFLAG) 

IF (ENDFLAG.EQ. ’END’ ) THEN 
WRITE (NTM. 2610) 

WRITE (NOT. 2610) 

RETURN 

60  ENDIF 

2610  FORMAT (’  *•*•  ERROR:  Insufficient  excitation  data;  at  least 
♦ two  values  must  be  given  for  excitation  functions.’) 

C 

C--7.0  CALL  PREDI  TO  DO  THE  WORK 
C 

CALL  PREDI (ID, K, C, EO, IA (MPT) .CH.NPR. IA (MPIEFN) . I A (MPEFN) . E. 

. IA (MPTD) . IA (MPTDD) , NPRT, IPRT, IWRT, NEFN. NNOD. NEQN . MBAN . ERR) 

RETURN 

END 

C 

ICOND 

SUBROUTINE  ICOND (ID. T, NNOD. NEQN . ERR) 

C-SUB: ICOND  - READS  AND  GENERATES  INITIAL  CONDITIONS 
C 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/  NTM. NTR. NIN. NOT. ND1 . ND2 . ND3 . ND4 
COMMON  /PARC/  FIN (12) . EXT (3) 

C 

C ICOND:  DATA  4 COMMON  STORAGE 

C 

REAL *8  T (NEQN) 

INTEGER  ID (NNOD).  IJX(3) 

CHARACTER  ENDFLAG* 3.  READ* 5 , FN(12)*1 
LOGICAL* 1 ERR 

SAVE  /IOLIST/, /PARC/ 

C 

C--1.0  CHECK  TO  SEE  IF  READ -FROM- D I SK  OPTION  IS  REQUESTED 
C 

CALL  FREE 

CALL  FREEH ( * ’, READ, 5,1) 

IF (READ. EQ. ’READ-’ ) THEN 
CALL  FREETY 

C GET  FILENAME  AND  CREATE  EXTENSION 

C — STORE  CURRENT  FILENAME  TEMPORARILY 
DO  10  1-1,12 
10  FN (I ) - FIN (I) 

CALL  FREEH ('O’ .FIN, 1. 12) 

C OPEN  FILE  <f 1 lenarae) . TMP , MOVE  TO  EOF.  BACKSPACE  ONE  RECORD 

CALL  FOPEN (ND1. ’TMP  ’.’OLD’) 

12  READ (ND1 , END-14) 

GO  TO  12 

14  BACKSPACE  ND1 

C READ  INITIAL  CONDITIONS  FROM  LAST  RECORD 

READ (ND1)  TIME, (T(N) ,N-1, NEQN) 

CALL  FCLOSE(NDl) 

WRITE (NOT,  2100)  (FIN (N) ,N«1,12) , (EXT (N) ,N«1,3) 

WRITE (NOT. 2200) 

WRITE (NOT. 2400)  (N . T (ABS (ID (N) ) ) . N-l . NNOD) 

C — RESTORE  CURRENT  FILENAME 
DO  16  1-1,12 
16  FIN (I)  - FN (I ) 

RETURN 

ENDIF 

C 

C— 2.0  WRITE  INITIAL  CONDITIONS  HEADER 
C 

WRITE  (NOT,  2200) 

C 

C— 3.0  INTERPRET  LINE  OF  DATA 
C 

30  CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH  ( ’ ' .ENDFLAG. 3.1) 

IF (ENDFLAG.EQ. ’END' ) GO  TO  40 
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C DO  IT 

TEMP  - 0.0 

CALL  FREER  ( 'T  ’.  TEMP . 1 ) 

CALLFREEIC  MJK(1)#3) 

IF (I JK (2) . EQ. 0)  IJK(2)-IJK(1) 

IF  (I  JK  (3)  . EQ.  0)  IJK(3)=1 
DO  34  N-IJK(l)  ,IJK(2)  ,IJK(3) 

IF(N.LE.O.OR.N.GT.NNOD)  THEN 

WRITE  (NTM.  2300)  N 

WRITE  (NOT.  2300)  N 

ERR-. TRUE. 

GO  TO  36 

END  IF 

34  T (ABS (ID (N) ) ) - TEMP 

C GET  NEXT  LINE  OF  DATA 

36  CALL  TREE 
GO  TO  30 
C 

C— 4.0  REPORT  INITIAL  CONDITIONS  IF  NO  ERROR  ENCOUNTERED 
C 

40  IF (ERR)  RETURN 

WRITE  (NOT,  2400)  (N.  T (ABS  (ID  (N)  ) ) ,N=l,NNOD) 

RETURN 

2100  FORMAT (/'  — NOTE:  Initial  conditions  read  from  last  record 

of:  ' 


. 1 6A1 ) 

2200  FORMAT (/ 


. 6X, 'Node 

Temp.  Node  Temp . 

Node 

Temp 

. Node 

. Temp . 

2300  FORMAT (' 

Node  Temp  . ' ) 

••*•  ERROR:  (Generated) 

Node  ' 

. 15.  ' 

Is  out  of  range 

2400  FORMAT ( (6X, 5 (14 ,G1 1 . 3)  ) ) 

END 

c PREDI 

SUBROUTINE  PREDI  (ID, K. C, EO. T, CH . NPR . IEFN . EFN . E. TD. TDD, 

♦NPRT, IPRT. IWRT.NEFN.NNOD.NEQN.MBAN, ERR) 

C-SUB: PREDI  - THE  KERNEL  OF  PREDIC 

IMPLICIT  REAL*  8 (A-H.O-Z) 

C 

C CAL-SAP:  DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/NTM, NTR, NIN. NOT. ND1 , ND2 . ND3 . ND4 
C 

C L1STP : DATA  4 COMMON  STORAGE 

C 

REAL* 8 K (NEON, MBAN ) , C (NEQN . MBAN ) .EO(NNOD) , T (NEQN ) , CH  (NNOD)  . 
.EFN  (NEFN.2)  ,E(NEQN)  . TD  (NEQN)  .TDD  (NEQN) 

INTEGER  ID (NNOD),  NPR  (NPRT) , IEFN (NNOD. 2) 

LOGICAL*l  ERR.  SKIP,  TDOF 

COMMON  /LIST/  TIME (3 ), ALPHA. TOLD, TNEW 
SAVE  /IOLIST/, /LIST/ 

WRITE (NOT. 2000) 

WRITE (NTM. 2000) 

2000  FORMAT ( / ' — TIME  STEP*) 

C 

C— 1.0  COMPUTE  INITIAL  TEMPERATURE  RATES:  (dT(0)/dt) 

C 

C 1.1  FOR  TEMP -DOF:  SCALE  (C]=[C)*1E15  OR  SET  (C]*1E15  FOR 

C MASSLESS  NODES 

C 

SKIP  = .FALSE. 

DO  10  N-l.NEQN 
IF (TDOF  (N. ID.NNOO) ) THEN 
C (N , 1 ) - C (N,  1)  *1 . 0E1S 

IF (C (N, 1 ) . EQ .0.0)  C (N , 1 ) - 1.0E15 
ENDIF 

10  CONTINUE 
C 

C CHECK  FOR  STEP  CHANGE  IN  SPECIFIED  INITIAL  TEMP 

C 

DO  12  11*1,  NNOD 
NN  - ID  (N) 

IFOm.LT.O)  THEN 
WNN  * ABS  (NN) 

NFNT  - IEFN(N.l) 

TEMP  - 0.0 

I F (NFNT. WE. O)  TEMP  ■ EFN  (NFNT.  1) 

IF(T(NNN) .NE.TEMP)  THEN 

WRITE (NTM. 2100) 

WRITE (NOT,  2100) 

SKIP  - .TRUE. 

ENDIF 

ENDIF 

12  CONTINUE 

2100  FORMAT ( / ' — NOTE:  Unable  to  estimate  time  step  for  initial 

♦ step  change  In  specified  nodal  temperature (s) . ' ) 


C 

C 1.2  FORM  (E(t-0) ) 

C 

TM  - 0.0 

CALL  EXCIT (ID, C, T, TD, E. CH, EO, IEFN . EFN . TM. NPRT. IPRT.  IWRT, 

+NEFN . NN OO. NEQN . MBAN . ERR) 

IF (ERR)  RETURN 

C 

c 1.3  FORM  RHS  ( E) - ( K) (T)  FOR  FLUX-DOF,  ( dT/dt ) *DIAG [C)  FOR  TEMP -DOF 

C 

CALL  RHS (ID. T,TD. E. K.C.NNOO, NEQN, MBAN) 

C 

C 1.4  SOLVE  (C)  (dT/dt)  *=  (E) 

C 

CALL  SYMBC (C, TD, NEQN , MBAN , NEQN, 1 ) 

IF (SKIP)  GO  TO  30 
C 

C— 2.0  COMPUTE  TIMESTEP  CHECK 
C 

C 2.1  COMPUTE  INITIAL  RATE  OF  TEMP  RATES:  d/dt (dT  (0) /dt ) 

C 

CALL  RHS (ID. TD, TDD, TDD. K.C. NNOD, NEQN. MBAN) 

CALL  SYMBC (C. TDD. NEQN. MBAN. NEQN, 2) 

C 

C— 2.2  COMPUTE  NORMS:  ||(T(0))||.  | | (dT  (0) /dt ) | | , | I d/dt  (dT  (0)  /dt ) I I 

C 

TN  - 0.0 
TDN  - 0.0 
TDDN  =0.0 
DO  22  N-l.NEQN 
TN  - TN  ♦ T (N)  **2 
TDN  * TDN  ♦ TD (N) **2 
22  TDDN  « TDDN  ♦ TDD(N)**2 
TN  - SQRT(TN) 

TDN  - SQRT (TDN) 

TDDN  - SQRT (TDDN) 

C 

C 2.3  EVALUATE  TAYLORS  EXPRESSION  FOR  TIME  STEP  ESTIMATE 

C 

B - 0.05 

IF (TDDN. NE. 0.0)  THEN 

DTEST  = (B*TDN  + SQRT (B*B*TDN*TDN  ♦ 2 . 0*B*TN*TDDN) ) /TDDN 
WRITE (NTM, 2200)  B*100.0.  DTEST, TIME (3 ) 

WRITE (NOT, 2200)  B*100.0.  DTEST. TIME (3) 

2200  FORMAT (/'  — NOTE:  Estimated  time  step  to  limit  error  to 

. approx. ' ,F5. 2. '%  is:'.G10.3,/ 

. * Specified  time  step  ls:',G10.3) 

ELSE 

WRITE (NTM. 2210) 

WRITE (NOT, 2210) 

2210  FORMAT (/'  — NOTE:  Unable  to  estimate  time  step  to  limit 

. error  for  the  given  system.') 

ENDIF 

C 

C--3.0  READ  (C)  4 (K)  FROM  DISK 
C 

30  REWIND  ND1 

READ (ND1 ) ( (K  (N , M) . N=1 , NEQN) , M=1 , MBAN) . 

. ((C(N,M) ,N=1 .NEQN) , M=1 , MBAN) 

C 

C--4.0  FORM  LHS:  ((C)  ♦ aDT(K)]  SCALING  TEMP -DOF  DIAG  TERMS 
C 

ADT  = ALPHA* TIME (3) 

DTA  - TIME (3)  - ADT 
DO  40  N-l, NEQN 
DO  40  M-l.MBAN 

40  C(N.M)  * C(N.M)  ♦ ADT*K (N, M) 

DO  42  N-l.NEQN 

42  IF (TDOF (N. ID.NNOO))  C(N.l)  * C (N, 1) *1 . 0E15 
C 

C— 5.0  TIME  STEP  THRU  SOLUTION 
C 

WRITE (NOT.  2500) 

WRITE (NTM,  2500) 

2500  FORMAT (/'  — RESPONSE') 

IOP  - 1 
ISTEP  - 0 

DO  500  TM*TIME (1 ) +TIME (3) ,TIME(2) .TIME (3) 

ISTEP  - ISTEP  ♦ 1 
C 

C 5.1  FORM  (E) 

C 

CALL  EXCIT (ID, C, T, TD, E. CH, EO, IEFN . EFN . TM, NPRT. IPRT, IWRT. 

♦NEFN . NNOD . NEQN . MBAN . ERR ) 

IF (ERR)  RETURN 

C 

C 5.1  PARTIAL  UPDATE  OF  T:  (T)  - (T)  ♦ (1-a) DT( dT/dt ) 

C 

DO  51  N=1 . NEQN 
51  T (N)  - T(N)  ♦ DTA*TD  (N) 

C 

C 5.3  FORM  RHS: (E)-(K) (T)  FOR  FLUX-DOF,  (dT/dt ) *DIAG (C)  FOR  TEMP-DOF 

C 

CALL  RHS (ID. T.TD.E. K.C. NNOD. NEQN, MBAN) 

C 

C 5.4  SOLVE  FOR  (dT/dt) 

C 

CALL  SYMBC (C. TD. NEQN. MBAN, NEQN. IOP) 

IOP  * 2 
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c 

C 5.5  COMPLETE  UPDATE  OF  T:  (T)  * (T)  ♦ aDT{dT/dt) 

C 

DO  55  N*1.NEQN 
55  T(N)  - T (N)  + ADT*TD(N) 

C 

C 5.6  REPORT  RESULTS 

C 

IF (MOD (ISTEP,  IPRT)  . EQ. 0)  THEN 

C WRITE  TO  <FILENAME> . TMP  OPTION 

IF(IWRT.EQ.l)  THEN 

CALL  FOPEN  (ND1. ’TMP  * .’NEW’) 

WRITE (ND1 ) TM, (T (N) ,N=1 , NEQN) 

CALL  FCLOSE (ND1 ) 

END  IF 

C PRINTABLE  OUTPUT 

WRITE  (NOT,  2510)  TM 

WRITE  (NOT,  2520) 

WRITE  (NOT,  2530)  (NPR  (N)  . T (ABS  (ID  (NPR  (N) ) ) ) ,N*1  ,NPRT) 

WRITE  (NTM,  2510)  TM 

WRITE  (NTM.  2520) 

WRITE  (NTM.  2530)  (NPR  (N)  , T (ABS  (ID  (NPR  (N) ) ) ) ,N*=1  ,NPRT) 

END  IF 

500  CONTINUE 
RETURN 

2510  FORMAT (/ ' Time:  '.Gil. 3) 

2520  FORMAT ( 

.6X,'Node  Ten*>.  Node  Temp . Node  Ten*).  Node 
. Ten*>.  Node  Ten*>.’) 

2530  FORMAT (( 6X, 5(14, G1 1.3))) 

END 

c ECNTRL 

SUBROUTINE  ECNTRL (ID. IEFN. NNOO, ERR) 

C-SUB:  ECNTRL  - FORMS  EXCITATION  FUNCTION  NUMBER  ARRAY  IEFN  (NNOO. 2) 

C 

C CAL -SAP : DATA  4 COMMON  STORAGE 

C 

COMMON  /IOLIST/NTM. NTR . NIN, NOT. ND1 . ND2 , ND3 . ND4 
C 

C ECNTRL:  DATA  4 COMMON  STORAGE 

C 

INTEGER  ID(NNOD) . IEFN (NNOO. 2) . IJK(3) 

LOGICAL *1  ERR 
CHARACTER  END FLAG *3 

SAVE  /IOLIST/ 


CHARACTER  ENDFLAG*  3 
C 

C— 1.0  UPDATE  OLD  VALUES 
C 

TOLD  * THEN 
DO  10  I -1, NEFN 
10  EFN(I.l)  - EFN (1,2) 

C 

C--2.0  READ  NEW  VALUES 
C 

CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END" 

CALL  FREEH ('  ’. ENDFLAG . 3 . 1 ) 

IF (ENDFLAG. EQ. 'END' ) RETURN 
CALL  FREER  CT’ .TNEW,  1) 

CALL  FREER ( 'N ' , EFN (1,2) , NEFN) 

RETURN 

END 


C RHS 

SUBROUTINE  RHS (ID. T.  TD.  E. K. C. NNOD. NEQN. MBAN) 

C-SUB : RHS  - FORMS  RHS  OF  (C](dT/dt)  - (E*)  - (E)  - (K]{T); 

C 

C {EMt))  * [E(t)  ) - [K)  (T(t)  ) FOR  FLUX-DOF 

C (EMt))  - (dT(t) /dt)*DIAG  OF  [Cl  FOR  TEMP-DOF 

C 
C 

C <E*)  IS  WRITTEN  OVER  (TD) 


IMPLICIT  REAL* 8 (A-H.O-Z) 

REAL*  8 T (NEQN)  . TD  (NEQN)  , E (NEQN)  , K (NEQN . MBAN ) , C (NEQN . MBAN  ) 
INTEGER  ID (NNOD) 

LOG I CAL *1  TDOF 

DO  20  N=l, NEQN 

C SCALE  BY  DIAGONAL  FOR  TEMP  PRESCRIBED  NODES 

IF (TDOF (N. ID. NNOD))  THEN 
TD  (N)  * TD  (N)  *C(N,  1) 

C FORM  [E) -IK) (T ) WHERE  (K)  IS  IN  COMPACT  STORAGE 

ELSE 

TEMP  « E(N)  - K (N,  1)  *T  (N) 

KK*=N 

LL=N 

DO  10  1*2 . MBAN 


WRITE (NOT, 2000) 

2000  FORMAT (/*  — EXCITATION  FUNCTION  NUMBERS’// 

.6X'Node  Temp  Ext. Flux  Ext . Temp .' ) 

10  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END- 

CALL  FREEH ( 1 ’ .ENDFLAG. 3. 1) 

IF (ENDFLAG. EQ. ’END’ ) RETURN 

C INTERPRET  LINE  AND  GENERATE  DATA 

IQ  * 0 

CALL  FREEI  CQ'.XQ.l) 

IT  « 0 

CALL  FREEI ( 'T' , IT, 1) 

ITEX  = 0 

CALL  FREEI ('X’ , ITEX. 1) 

CALL  FREEI ('  * ,IJK(1),3) 

IF  (I  JK  (2)  . EQ.  0)  I JK  (2)  *IJK  (1) 

IF  (I  JK  (3)  . EQ.  0)  I JK  (3)  *1 
DO  12  N=I  JK  (1 ) , I JK  (2),IJK(3) 

IF (N.LE.O.OR.N.GT.NNOD)  THEN 
WRITE  (NTM,  2010)  N 
WRITE  (NOT.  2010)  N 
ERR*. TRUE . 

GO  TO  10 

2010  FORMAT ( ' ****  ERROR:  (Generated)  Node’, 15,'  la  out  of  range.') 
ENDIF 
NN  - ID  (N) 

NNN  - ABS  (NN) 

C TEMPERATURE-PRESCRIBED  NOOES 

IF(NN.LT.O)  THEN 

WRITE  (NOT,  2020)  N.IT 


KK  * KK-1 

IF(KK.GT.O)  TEMP  = TEMP  - K (KK , I ) *T (KK) 

LL  * LL+1 

IF  (LL . LE  .NEQN)  TEMP  « TEMP  - K(N,I)*T(LL) 
10  CONTINUE 

TD  (N)  * TEMP 
ENDIF 

20  CONTINUE 
RETURN 
END 


C EXCIT 

SUBROUTINE  EXCIT (ID.C.T, TD, E.CH.EO, IEFN . EFN . TM. NPRT . IPRT. IWRT, 
♦NEFN , NNOD, NEQN . MBAN . ERR ) 

C--SUB : EXCIT  - FOR  t-TIME; 

C FOR  FLUX-DOF  COMPUTES  (E(t)); 

C 

C (E  (t ) > - (EO)  ♦ (Q(t)  ) ♦ (CH)T(TEX(t)  ) 

C 

C FOR  TEMP-DOF  COMPUTES  (dT/dt(t)} 

C 

C (dT(t)/dt)  - (l/Dt)(T(t)  - T(t-Dt)) 

C 

C FROM  EXCITATION  FUNCTION  DATA  STORED  AND  READ  INTO  EFN (NEFN. 2) 

C AND  LIST  OF  FUNCTION  NUMBERS 


IMPLICIT  REAL* 8 (A-H.O-Z) 


IEFN(N.l)  - IT 
2020  FORMAT (CX, 14. 10X, 15) 

C FLUX-PRESCRIBED  NOOES 

ELSEIF(NN.GT.O)  THEN 

WRITE (NOT. 2030)  M, IQ. ITEX 
IEFN  (N, 1 ) - IQ 

IEFN (N, 2)  - ITEX 
2030  FORMAT ( 6X , 1 4 , 2 OX, 1 5 , 1 OX , 1 5 ) 
ENDIF 

12  CONTINUE 
GO  TO  10 


C 

c 

c 

c 

c 

c 


CAL -SAP : DATA  4 COMMON  STORAGE 

COMMON  /IOLIST/NTM, NTR, NIN. NOT, ND1, ND2, ND3, ND4 

EXCIT:  DATA  6 COWK3N  STORAGE 

REAL* 0 C (NEQN , MBAN) . T (NEQN) , TD (NEQN) . CH (NNOD) , EO  (NNOD) , E (NEQN) , 
♦EFN (NEFN, 2) 

INTEGER  ID (NNOO).  IEFN  (NNOO. 2) 

LOGICAL* 1 ERR 
CHARACTER  ENDFLAG* 3 


END 

C GETEFN 

SUBROUTINE  GETEFN (EFN. TOLD, TNEW, NEFN, ENDFLAG) 

C-SUB :GETEFN  - READ  INPUT  FILE  TO  UPDATE  EXCITATION  FUNCTION  DATA 

IMPLICIT  REAL* 8 (A-H.O-Z) 


COMMON  /LIST/  TIME (3) .ALPHA, TOLD, TNEW 
SAVE  /IOLIST/, /LIST/ 


C 

C--1.0  UPDATE  EXCITATION  FUNCTION  DATA  IF  NEEDED 
C 

10  IF (TM.GT.TNEW)  THEN 
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REAL* 8 EFN (NEFN. 2) 


CALL  GETEFN (EFN, TOLD, TNEW, NEFN, ENDFLAG) 


C--2.0  GET  PRINT  INTERVAL.  WRITE  OPTION 


C 

IF (ENDFLAG. EQ. 'END')  THEN 
WRITE (NTH. 2100) 

WRITE  (WOT,  2100) 

2100  FORMAT ( ' ••••  ERROR:  Insufficient  excitation  data.') 

ERR  - .TRUE. 

RETURN 

ENDIF 
GO  TO  10 

ELSEIF (TM.LT.TOLD)  THEN 
WRITE  (NTH.  2110)  TM 

WRITE  (NOT,  2110)  TM 

2110  FORMAT (•  *•**  ERROR:  Excitation  data  not  defined  for  current 
. time: * .G10.3) 

ERR-. TRUE. 

RETURN 

ENDIF 

C INTERPOLATION  FRACTION 

XT  ■ (TM-TOLD) / (TNEW-TOLD) 

C 

C— 2.0  FORM  (E)  OR  (dT/dt) 

C 

CALL  SEROR (E. NEQN . 1 ) 

DO  20  N-l.NNOD 
NN  « ID  (N ) 

NNN  - ABS  (NN) 

C FLUX-DOF:  SUM  INT.  FLUX  + EXT.  FLUX  + EFFECT.  CONVECTIVE  FLUX 

IF(NN.GT.O)  THEN 
NFNQ  = I EFN  (NNN. 1) 

NFNC  * I EFN (NNN. 2) 

Q - 0.0 

IF  (NFNQ.NE.  0)  Q - XT*  (EFN  (NFNQ,  2) -EFN  (NFNQ,  1 ) ) 4-  EFN  (NFNQ.  1) 

CONVE  = 0.0 

IF  (NFNC.NE.  0)  CONVE  - CH  (N)  • (XT*  (EFN  (NFNC . 2)  - EFN  (NFNC.  1)) 

4-  ♦ EFN  (NFNC.  1)) 

E (NNN)  - E(NNN)  + EO(N)  ♦ Q + CONVE 

C TEMP -DOF:  COMPUTE  (dT/dt)  FROM  PRESCRIBED  TEMP 

ELSEIF (NN.LT. 0)  TEEN 
MFNT  * IEFN(NNN.l) 

TEMP  -0.0 

IF(NFNT.NE.O)  TEMP  «=  XT*  (EFN  (NFNT,  2) -EFN  (NFNT,  1)  ) 4-  EFN(NFNT.l) 

TD (NNN)  - (TEMP  - T (NNN) ) /TIME (3) 

ENDIF 

20  CONTINUE 

RETURN 
END 

C REPORT 

SUBROUTINE  REPORT  (MPNPR,  NPRT,  IPRT,  IWRT.  NNOD,  ERR) 

C-SUB:  REPORT  - READS  REPORT  CONTROL  INFORMATION 
C FORMS  REPORT  CONTROL  ARRAY  NPR  (NPRT) 

C 

e CAL -SAP : DATA  & COMMON  STORAGE 

C 

CHARACTER  EXT*1,  FIN*1 
COMMON  MTOT.NP, IA(1) 

COMMON  /DBSYS/  NUMA.NEXT, IDIR, IP (3) 

COMMON  /IOLIST/  NTM, NTR, NIN . NOT, ND1 . ND2 , ND3 , WD4 
COMMON  /PARC/  FIN (12) , EXT (80) 

SAVE  /DBSYS/. /IOLIST/, /PARC/ 

C 

C- REPORT:  DATA  4 COMMON  STORAGE 

C 

INTEGER  I JK (3) 

LOGICAL* 1 ERR 

CHARACTER  ENDFLAG* 3,  WRITEFLAG*5 

SAVE  /IOLIST/ 

C 

C— 0.0  WRITE  EEAfiER 
C 

WRITE (NTM, 2000) 

2000  FORMAT  (/’  — REPORT  CONTROL’) 

C 

C— 1.0  FIND  "REPORT" 

C 

CALL  FINDN  ( ’REPORT’  , 6,  KEY) 

IF (KEY . EQ. 1)  THEN 
WRITE  (NTM.  2100) 

WRITE  (NOT.  2100) 

ERR= . TRUE . 

RETURN 
ENDIF 

CALL  FREETY 
C 


CALL  FREE 
CALL  FREETY 

CALL  FREEH (’  * , WRITEFLAG. 5 , 1 ) 

IWRT— 0 

IF (WRITEFLAG. EQ. ’WRITE’ ) THEN 
IWRT  - 1 
WRITE  (NTM, 2200) 

WRITE (NOT. 2200) 

ENDIF 

IPRT-1 

CALL  FREEI ( ’ T’ , IPRT,  1 ) 

C 

C--3.0  FORM  PRINT  CONTROL  ARRAY 
C 

C DEFINE  ARRAY  NPR  (1.1)  TO  GET  POINTER  4 CREATE  DIRECTORY  ENTRY 

CALL  DELETE ( ’NPR  ’) 

CALL  DEFINE (’NPR  ’.MPNPR.1,1) 

C GENERATE  PRINT  CONTROL  ARRAY  WHILE  BLANK  COMMON  SPACE  AVAILABLE 

NPRT— 0 

30  CALL  FREE 
CALL  FREETY 

C CHECK  FOR  "END" 

CALL  FREEH (’  ’, ENDFLAG, 3 . 1 ) 

IF (ENDFLAG. EQ. ’END’ ) GO  TO  36 
C DO  IT 

CALLFREEK’  •,XJK<1),3) 

IF  (I  JK  (2)  .EQ.  0)  I JK  (2)  -I  JK  (1 ) 

IF  (I  JK  (3)  . EQ.  0)  I JK  (3)  -1 
DO  34  N-IJK  (1)  , I JK  (2) , I JK  (3) 

IF(N.LE.O.OR.N.GT.NNOD)  THEN 
WRITE (NTM. 2300)  N 
WRITE (NOT. 2300)  N 
ERR-. TRUE. 

GO  TO  30 
ENDIF 

NN  - HPNPR+N-1 

C CHECK  BALNK  COMMON  SPACE 

IF (NN . GT . IDIR)  THEN 
WRITE (NTM, 2310) 

WRITE (NOT. 2310) 

ERR= . TRUE . 

RETURN 

ENDIF 

NPRT  - NPRT  4-  1 
3 4 IA  (NN)  - N 
GO  TO  30 

C REDEFINE  DIRECTORY  FOR  ACTUAL  SIZE  OF  NPR  ARRAY 

36  CALL  DELETE  (’NPR  ’) 

CALL  DEFINE ( ’ NPR  ’ . MPNPR , NPRT, 1 ) 

RETURN 

2100  FORMAT ( ’ •***  ERROR : No  report  control  Information  found.’) 

2200  FORMAT ( 

.’  — Solution  will  be  written  to  binary  file  <f llename> . TMP 

. In  form  (TIME. T (N) ,N=l,NNOO) . ’ ) 

2300  FORMAT ( ’ ••*•  ERROR:  (Generated)  Node, ’15*  selected  for 
printable 

. output  Is  out  of  range.’)  , 

2310  FORMAT ( » 

.’  ••**  ERROR:  Insufficient  storage  to  form  print  control  array. 
.*  Increase  IA (MTOT)  by  ’,16) 

END 

C 

TDOF 

FUNCTION  TDOF (NEQ, ID. NNOD) 

C-FUN : TDOF  - DETERMINES  IF  EQUATION  NUMBER  NEQ  IS  A TEMP  DOF 
LOGICAL* 1 TDOF 
INTEGER  ID (NNOD) 

TDOF  « .FALSE. 

DO  10  N-l.NNOO 

IF  ( (ID  (N) . LT . 0)  .AND.  (ABS (ID (N) ) . EQ. NEQ) ) THEN 
TDOF  - .TRUE. 

RETURN 

ENDIF 

10  CONTINUE 
RETURN 
END 

C* ••••*•• *••••••• •*••»**•* ft ft**************** ft •••••••••* • •• ******** 

C CAL  - SAP  LIBRARY  EXTENSIONS 
(;**••*•****•***••••*••••*••**»••*••***••**•••••••••*•*•»•*•*•••**••** 

c FREETY 

SUBROUTINE  FREETY 

C— SUB:FREETY  - TYPE  COMMAND  LINE  TO  SCREEN 
CHARACTER* 1 LINE 
COMMON  /ILINE/  II 
COMMON  /CLINE/  LINE (160) 

COMMON  /IOLIST/  NTM, NTR, NIN . NOT. ND1 . ND2 . ND3 , ND4 
SAVE  /ILINE/. /CLINE/. /IOLIST/ 

WRITE (NTM. 2000)  (LINE (I) , 1*1. II) 

RETURN 

2000  FORMAT  (1X.8QA1) 

END 
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c FINDN 

SUBROUTINE  FINDN (SEP, NC, KEY) 

CHARACTER* 1 LINE. SEP (NC) 

COMMON  /CLINE/  LINE  (160) 

COMMON  /ILINE/  II 

COMMON  /IOLIST/  NTM, NTR, NIN , NOT, ND1 , ND2 , ND3 , ND4 
SAVE  /CLINE/. /ILINE/, /IOLIST/ 

C FIND  SEPARATOR  OF  NC  CHARACTERS  IN  INPUT  FILE 

KEY=0 

REWIND  NIN 
50  CONTINUE 

READ (NIN, 1000, ERR=200, END=200)  (LINE(I) ,I=1,NC) 

II  * NC 
CALL  UPPER 
DO  60  N-l.NC 

60  IF (SEP (N) .NE. LINE (N) ) GO  TO  50 
GO  TO  900 

200  KEY=1 
900  RETURN 

1000  FORMAT  (80A11 
END 

C SAVE 

SUBROUTINE  SAVE 

C— SUB:  SAVE  - WRITES  INCORE  DATA  BASE  TO  <FILENAME> . DBS 
CHARACTER  EXT*1.  FIN*1 
COMMON  MTOT. NP.,  IA  (10000) 

COMMON  /DBSYS/  NUMA, NEXT. IDIR, IP (3) 

COMMON  /IOLIST/  NTM. NTR. NIN . NOT. ND1 , ND2 , ND3 , ND4 
COMMON  /PARC/  FIN (12) . EXT (80) 

SAVE  /DBSYS/, /IOLIST/. /PARC/ 

CALL  FOPEN (ND1, 'DBS  • .'NEW') 

WRITE (ND1)  NUMA. NEXT,  IDIR.  MTOT. IP 

WRITE (ND1)  (IA (I) , 1=1 , NEXT) . (IA ( J)  . J=IDIR,  MTOT) 

CALL  FCLOSECNDl) 

RETURN 

END 

c ZEROI 

SUBROUTINE  ZEROI (IA. NR. NC) 

DIMENSION  IA (NR, NC) 

DO  10  1=1, NR 
DO  10  J=1 , NC 
IA  (I,  J)  = 0 
10  CONTINUE 
RETURN 
END 

C ZEROR 

SUBROUTINE  ZEROR (A. NR. NC) 

REAL*  8 A (NR,NC) 

DO  10  1=1, NR 
DO  10  J=1 , NC 
A (I , J)  = 0.0 
10  CONTINUE 
RETURN 
END 

c SYMBC 

SUBROUTINE  SYMBC (C . B. MAX , MB . ND . IOP ) 

C— SUB:  SYMBC  - SYWETRIC  COMPACT-FORM  EQUATION  SOLVER 
C SOLVES  SYMMETRIC  BANDED  EQUATIONS; 

C 

C [C)  (X)  = ( B) 

C 

C IN-CORE  WITH  [C]  STORED  IN  COMPACT  FORM  C (ND, MB) . 

C SOLUTION  X (ND)  IS  WRITTEN  OVER  B (ND)  . 

C 

C MAX  = MAXIMUM  DOF  FOR  PARTIAL  L*D*U  DECOMPOSITION 

C MB  = (HALF)  BANDWIDTH 

C ND  = NUMBER  OF  DEGREES  OF  FREEDOM 

C IOP  = 1 : COMPLETE  SOLVE 

C IOP  * 2 : RESOLVE 

C 


IMPLICIT  REAL*  8 (A-H.O-Z) 

DIMENSION  C (ND,  MB)  , B (ND) 

C-l — PERFORM  THE  L*D*U  DECOMPOSITION  OF  C AND  FORM  Y 
NBM  - MB- 1 
DO  300  M*2 , MAX 

N - M-l 

ir(C(N,l) .EQ.0.0)  OO  TO  300 
WN  - B(N) 

NNB  - MIN  ( (N+NBM)  .MAX) 

KJ  - NNB  - N + 1 
II  - 1 

DO  200  I - M, NNB 

XX  ■ XX  ♦ 1 

IF  (C (N, II) .EQ.0.0)  GO  TO  200 
GIN  - C (N,  II)  /C  (N,  1) 

B (I)  = B (I ) - GIN*WN 
IF (IOP . EQ. 2)  GO  TO  200 
J = 0 

DO  100  K-II.KJ 
J = J + 1 

100  C(I,J)  = C(I,J)  - C (N , K)  *GIN 
200  CONTINUE 
300  CONTINUE 

C -2 --PERFORM  BACKSUBSTITUTION 


M » MAX 

400  IF (C  (M, 1) .NE.0.0)  B (M)  «B(M)/C(M.l) 

M = M - 1 

IF(M.LE.O)  GO  TO  600 
KJ  = MIN  (MB. MAX-M+1) 

J = M 

DO  500  N=2 , KJ 
J = J + 1 

500  B(M)  = B (M)  - C (M,  N ) *B  ( J) 

GO  TO  400 
600  RETURN 
END 

C SYMBCX 

SUBROUTINE  SYMBCX (C. B. MAX. MB. ND, IOP) 

C — SUB: SYMBCX  - SY METRIC  COMPLEX  COMPACT-FORM  EQUATION  SOLVER 
C SOLVES  SYMMETRIC  BANDED  COMPLEX  EQUATIONS; 

C 

C (C*J  (X*)  = (B*) 

c 

C IN-CORE  WITH  (C*)  STORED  IN  COMPACT  FORM  C (ND. MB) . 

C SOLUTION  X*  (ND)  IS  WRITTEN  OVER  B* (ND)  . 

C 

C MAX  » MAXIMUM  DOF  FOR  PARTIAL  L*D*U  DECOMPOSITION 

C MB  - (HALF)  BANDWIDTH 

C ND  » NUMBER  OF  DEGREES  OF  FREEDOM 

C IOP  - 1 : COMPLETE  SOLVE 

C IOP  - 2 : RESOLVE 

C 

IMPLICIT  COMPLEX  (A-H.O-Z) 

DIMENSION  C (ND. MB) , B (ND) 

C-l— PERFORM  THE  L*D*U  DECOMPOSITION  OF  C AND  FORM  Y 
NBM  = MB - 1 
DO  300  M=2 , MAX 
N = M-l 

IF  (C  (N , 1 ) . EQ . ( 0 . 0 , 0 . 0 ) ) GO  TO  300 
WN  = B (N) 

NNB  = MIN( (N+NBM) .MAX) 

KJ  = NNB  - N + 1 
II  = 1 

DO  200  I = M. NNB 
XI  - XX  + 1 

IF(C(N.II)  .EQ.  (0.0.  0.0)  ) GO  TO  200 
GIN  = C(N.II)/C(N,1) 

B (I ) = B (I ) - GIN*WN 
IF (IOP . EQ. 2 ) GO  TO  200 
J = 0 

DO  100  K=II , KJ 
J « J + 1 

100  C(I.J)  = C(I,J)  - C (N,  K)  *GIN 
200  CONTINUE 
300  CONTINUE 

C-2- -PERFORM  BACKSUBSTITUTION 
M = MAX 

400  IF(C(M.  1)  .NE.  (0.0,0. 0)  ) B (M)  =B(M)/C(M,1) 

M = M - 1 

IF(M.LE.O)  GO  TO  600 
KJ  = MIN (MB. MAX-M+1) 

J - M 

DO  500  N-2.KJ 
J = J + 1 

500  B (M)  = B (M)  - C (M, N)  *B  ( J) 

GO  TO  400 
600  RETURN 
END 

CROUT 

SUBROUTINE  CROUT (A, B. N, ND, LD. KEY) 

C 

C GENERAL  EQUATION  SOLVER  - M.I.  HOIT  1/17/83 
C SYMMETRIC  AND  NON- SYMMETRIC 

C 

C KEY  - 0 TRIANGULARIZE  AND  SOLVE 

C KEY  = 1 TRIANGULARIZE  ONLY 

C KEY  - 2 FORWARD  AND  BACK  SUBSTITUTION  ONLY 

C 

IMPLICIT  REAL* 8 (A-H.  O-Z) 

DIMENSION  A (ND, ND) , B (ND, LD) 

C 

II  - 1 

IF (A(l.l) .NE.0.0)  GO  TO  10 

PAUSE  ' •***  ERROR: SUB: CROUT:  Attenpt  to  solve  singular  system. 
GO  TO  900 
10  IM  = 1 

IF (KEY. EQ. 2)  GO  TO  200 
IF (N . EQ . 1 ) GO  TO  200 
DO  140  J-2.N 
JM  - J - 1 

A ( J,  1)  - A ( J,  1 ) /A  (1 , 1)  “ 

IF ( J .EQ. 2)  GO  TO  130 
DO  120  1-2. JM 
IM  - I - 1 
DO  100  K=1 , IM 

C FORM  L (I,  J)  4 U(I.J) 

A(J.I)  = A ( J,  I)  - A ( J,  K)  *A  (K,  I) 

100  A(I.J)  = A (I , J)  - A (I , K)  *A  (K.  J) 

D = A (I , I) 

C CHECK  IF  A ZERO  IS  ON  DIAGONAL 

IF(D. NE.0.0)  GO  TO  120 

PAUSE  * *•**  ERROR: SUB: CROUT:  Attempt  to  solve  singular  system. 
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GO  TO  §00 

120  A(J.I>  - A(J,I>  /D 
130  DO  140  K-1,JH 

140  A(J.J)  - A(J.J)  - A ( J,  K)  • A <K,  J) 

C FORWARD  AND  BACK  SUBSTITUTE 

200  IF (KEY . EO . 1 ) GO  TO  900 
DO  400  L-l.LD 

C FORM  T (1*  L) 

IF(N.EQ.l)  GO  TO  300 
DO  210  1-2. N 
IN  - I - 1 
DO  210  K-l.IH 

210  B(I,L>  - B(I.L)  - A (I , K)  *B  (K,  L) 

C FORM  X (X«  L) 

300  B(N,L)  - B (N.L)  /A(N.N) 

IFOI.EO.l)  GO  TO  400 

IN  - * - 1 

DO  330  1-2. N 
JH  « IN  ♦ 1 
DO  310  J « JM.N 

310  B (IN,  L)  - B(IH.L)  - A (IM.  J)  *B  ( J.  L) 
320  B(IH.L)  - B (IN.  L) /A  (IN.  IN) 

330  IN  - IN  - 1 
C 

400  CONTINUE 
900  RETURN 

END 
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Building  Energy  Simulation  that  are  expected  to  provide  a means  to  unify  existing 
building  energy  simulation  theory.  DTAM1  provides  a library  of  discrete  thermal 
elements,  that  may  be  assembled  to  model  thermal  systems  idealized  to  have  constant 
material  and  heat  transfer  properties  (i.e.,  linear  idealizations),  including: 

- ID  two-node  thermal  resistance  elements 

- single-node  lumped  capacitance  elements 

- two-node  fluid  flow  loop  element 

- ID  two-to-four  node  isoparametric  conduction  Finite  Elements 

- 2D  four-node  isoparametric  conduction  Finite  Elements  (planar  and  axisymmetric) 
Equations  defining  a variable  node  mean  radiant  temperature  element  are  also  presented 
in  the  report.  Steady  state  and  transient  analysis  capabilities  are  included. 
Temperature,  heat  flow  rate,  and  convective  boundary  conditions  may  be  modeled  and 
s-stem  temperature  variables  may  be  constrained  to  be  equal  so  that  mixed  assemblages 
of  ID  and  2D  elements  may  be  employed. 
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